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Requirements  for  the  Degree  of  Doctor  of  Philosophy 

STUDY  OF  TMGa-NH3-N2  SYSTEM 
USING  in  situ  RAMAN  SPECTROSCOPY 

By 

Min  Huang 
May  2003 

Chair:  Dr.  Timothy  James  Anderson 
Major  Department:  Chemical  Engineering 

The  rapid  growth  of  the  semiconductor  industry  has  demanded  the  shortening  of 
process  development  time,  particularly  in  chemical  vapor  deposition  (CVD)  technology. 
This  poses  challenges  on  two  major  fronts:  developing  new  CVD  chemical  precursors 
and  optimizating  reactor  designs  and  growth  conditions.  The  objective  of  this 
dissertation  is  to  demonstrate  the  combination  of  in  situ  Raman  spectroscopy,  ab  initio 
calculations  for  product  properties,  numerical  CVD  reactor  modeling  and  inverse 
problem  solution,  to  identify  reaction  mechanisms,  reaction  rate  constants  and  species 
diffusivities. 

3-Dimentional,  spatially  resolved  gas  phase  temperature  and  species 
concentration  profiles  were  obtained  by  using  in  situ  Raman  spectroscopy  in  an  inverted, 
vertical  impinging  flow  reactor.  Analytical  methods  for  were  developed  to  reduce 
experimental  data.  For  concentration  measurement,  modification  of  the  "Relative 
normalized  differential  Raman  scattering  cross  section  “Zj"  formula  was  proposed. 
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An  inverse  problem  solution  code  was  developed  based  on  Tikhonov 
regularization  for  identifying  species  diffusivity  and  gas  phase  homogeneous  reaction  rate 
constants.  The  gas  phase  binary  mass  diffusion  coefficient  of  TMGa  in  N2  was  found  to 
be  0.08  [cm^s"'].  The  rate  constant  of  TMGa  homogeneous  decomposition  in  N2 
environment  was  found  to  be  3.4xl0'^  exp(60.1  kcal/RT)  [s"'].  The  rate  constant  of  NH3 
homogeneous  decomposition  in  N2  environment  was  found  to  be 
1.4xl0*^exp(90. 11  kcal/RT)  [cm^  mof’  s'*].  The  identified  values  matches  literature 
values  very  well. 

For  the  first  time,  TMGa-NH3-N2  reaction  system  was  studied  at  conditions  close 
to  the  actual  GaN  MOCVD  conditions.  A gas  phase  reaction  mechanism  was  proposed 
based  on  measured  gas  phase  temperature  and  the  species  concentration  profiles  of  the 
key  components  in  the  TMGa-NH3-N2  reaction  system. 

Given  the  proposed  TMGa-NH3-N2  reaction  mechanism,  a scale  analysis  and 
operating  window  calculation  for  GaN  MOCVD  was  performed  to  a provided  qualitative 
understanding  of  the  growth  process.  An  analytical  form  of  the  GaN  growth  rate  was 
deduced.  An  optimized  operating  window  for  GaN  growth  was  also  suggested.  Both 
agree  very  well  with  the  operating  conditions  provided  by  commercial  reactor 
manufacturers. 

Numerical  simulation  of  TMGa:NH3-NH3-N2  gas  phase  reaction  system  was 
successfully  carried  out.  Simulated  gas  phase  temperature  and  NH3  concentration 
profiles  agree  well  with  the  in  situ  Raman  experimental  results.  The  simulated 
concentration  profiles  of  the  reaction  products  were  in  agreement  with  of  the  in  situ 
Raman  experimental  results. 
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CHAPTER  1 
INTRODUCTION 

Philosophical  Background 

The  Definition  of  Truth 

The  ancient  Greek  philosopher  Aristotle  gave  his  definition  of  truth  in 
Metaphysics  [1]:  “To  say  of  what  is  that  it  is  not,  or  of  what  is  not  that  it  is,  is  false,  while 
to  say  of  what  is  that  it  is,  or  of  what  is  not  that  it  is  not,  is  true.”  Here  Aristotle 
suggested  that  the  truth  relies  on  the  observer’s  intuition.  Following  this  concept,  Alfred 
Tarski  presented  his  modem  day  definition  of  truth  [2]: 

“Snow  is  white”  is  tme  if,  and  only  if,  snow  is  white. 

Here  the  tmth  definition  itself  was  to  be  a definition  of  “True”  in  terms  of  the  other 
expressions  of  the  metalanguage.  So  the  definition  was  to  be  in  terms  of  syntax,  set 
theory  and  the  notions  expressible  in  L,  but  not  semantic  notions  like  “denote”  or 
“mean.” 

In  addition  to  intuition,  Tarski  pointed  out  that  the  definition  of  “Tme”  should 
also  be  presented  ‘formally  correct’  [2].  This  means  that  it  should  be  a sentence  of  the 
form: 

For  all  L,  Tme  if  and  only  if  Z,’,  L-1 

where  “Tme”  should  never  occur  in  L,  the  definition  should  be  provably  equivalent  to  a 
sentence  of  this  form.  The  equivalence  must  be  provable  using  axioms  of  the 
metalanguage  that  do  not  contain  “True.” 
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Tarski  also  pointed  out  that  the  definition  should  be  ‘materially  adequate’  [2], 

This  means  that  the  objects  satisfying  L’  should  be  exactly  the  objects  that  one  would 
intuitively  count  as  being  true  sentences  of  L,  and  that  this  fact  should  be  provable  from 
the  axioms  of  the  metalanguage  that  do  not  contain  “True.”  Metalanguage  is  one  in 
which  the  properties  of  a language  under  study,  the  object  language,  are  stated. 
Metalanguage  may  be  identical  to  the  object  language,  as  when  the  grammatical 
properties  of  English  are  stated  in  English,  but  it  is  often  distinct  from  the  object 
language. 

It  may  look  like  a paradoxical  requirement:  if  it  can  be  proven  what  Tarski  asked 
for  just  from  the  axioms  of  the  metalanguage,  then  it  must  be  already  have  a materially 
adequate  formalization  of  ‘true  sentence  of  L’  within  the  metalanguage.  Tarski  escaped 
that  by  using  (in  general)  infinitely  many  sentences  of  L’  to  express  truth  [2],  namely  all 
the  sentences  of  the  form: 

y/  (5)  if  and  only  if  ^ L-2 

whenever  5 is  the  name  of  a sentence  S of  L and  IPis  the  copy  of  S in  the  metalanguage. 
So  the  technical  problem  is  to  find  a single  formula  ^^that  allows  deducing  all  these 
sentences  from  the  axioms  of  L’.  This  formula  ^will  serve  to  give  the  explicit  definition 
of  True. 

Alfred  Tarski  proved  [2]  that: 

• Suppose  L=L’,  for  any  proposed  definition  of  truth,  one  of  the  adequacy  sentences 
can  be  outright  disproved  in  L.  Hence  no  adequate  definition  of  truth  in  L can  be 
given  in  L 

• Suppose  L’  includes  second-order  logic  over  L and  PA  (first  order  number  theory),  an 
adequate  definition  of  truth  in  L can  be  given  in  L’ 
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This  Alfred  Tarski’s  early  definition  of  truth  although  can  cause  the  so  called  “Liar 
Paradox”,  therefore,  problematic  in  some  senses  [3,4],  it  has  profound  significances  in 
shaping  methodologies  of  deductive  sciences.  The  definition  and  its  analogue  later 
directly  evolved  into  the  so  called  model  theory  and  the  definition  of  model-theoretic 
truth. 

Model  Theory 

To  model  a phenomenon  is  to  construct  a formal  theory  that  describes  and 
explains  it.  In  a closely  related  sense,  model  a system  or  structure  that  planed  to  build,  by 
writing  a description  of  it.  These  are  very  different  senses  of  ‘model’  from  that  in  the 
model  theory  here:  the  ‘model’  of  the  phenomenon  or  the  system  is  not  a structure  but  a 
theory. 

Model  theory  is  the  study  of  the  interpretation  of  any  language,  formal  or  natural, 
by  means  of  set-theoretic  structures,  using  the  analogue  of  Alfred  Tarski's  truth 
definition.  Mainstream  model  theory  is  now  a sophisticated  branch  of  mathematics 
called  first  order  model  theory. 

Tarski’s  model-theoretic  truth.  When  a sentence  S being  written  or  spoken  in  the 
way  that  it  expresses  nothing  either  true  or  false,  because  some  crucial  information  is 
missing  about  what  the  words  mean.  If  this  information  was  added  on  so  that  S comes  to 
express  a true  or  false  statement,  it  may  be  said  to  interpret  S,  and  the  added  information 
is  called  an  interpretation  of  S.  If  the  interpretation  / happens  to  make  S state  something 
true,  it  may  be  said  that  / is  a model  of  S,  or  that  / satisfies  S,  in  symbols  ‘/  1=  5”.  Another 
way  of  saying  that  / is  a model  of  S is  to  say  that  S is  true  in  I,  and  the  notion  of  model- 
theoretic  truth  may  be  used,  which  is  truth  in  a particular  interpretation.  However,  the 
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statement 'S  is  true  in  F is  just  a paraphrase  of  ‘5,  when  interpreted  as  in  1,  is  true’;  so 
model-theoretic  truth  is  parasitic  on  plain  ordinary  truth.  Here  Tarski’s  model-theoretic 
truth  definition  is  biconditional,  necessary  and  sufficient,  and  therefore,  must  be  unique. 

This  interpretation  explains  (i)  what  objects  some  expressions  refer  to,  and  (ii) 
what  classes  some  quantifiers  range  over. 

Structure.  Interpretations  that  consist  of  items  (i)  and  (ii)  appear  very  often  in 
model  theory,  and  they  are  known  as  structures.  Particular  kinds  of  model  theory  use 
particular  kinds  of  structure;  for  example  mathematical  model  theory  tends  to  use  so- 
called  first-order  structures.  The  objects  and  classes  in  a structure  carry  labels  that  steer 
them  to  the  right  expressions  in  the  sentence.  These  labels  are  an  essential  part  of  the 
structure.  If  the  same  class  is  used  to  interpret  all  quantifiers,  the  class  is  called  the 
domain  or  universe  of  the  structure.  But  sometimes  there  are  quantifiers  ranging  over 
different  classes. 

Signature.  Each  mathematical  structure  is  tied  to  a particular  first-order  language. 
A structure  contains  interpretations  of  certain  predicate,  function  and  constant  symbols; 
each  predicate  or  function  symbol  has  a fixed  arity.  The  collection  K of  these  symbols  is 
called  the  signature  of  the  structure.  Symbols  in  the  signature  are  often  called  nonlogical 
constants.  The  first-order  language  of  signature  K is  the  first-order  language  built  up 
using  the  symbols  in  K,  together  with  the  equality  sign  =,  to  build  up  its  atomic  formulas. 
If  Ai  is  a signature,  S is  a sentence  of  the  language  of  signature  K and  ^4  is  a structure 
whose  signature  is  K,  then  because  the  symbols  match  up  then  A makes  S either  true  or 


false. 
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First  Order  Model  Theory 

First-order  model  theory,  also  known  as  classical  model  theory,  is  a branch  of 
mathematics  that  deals  with  the  relationships  between  descriptions  in  first-order 
languages  and  the  structures  that  satisfy  these  descriptions. 

The  elements  of^  are  the  elements  of  dom(^).  Likewise  the  cardinality  or  power 
of  A is  the  cardinality  of  its  domain.  If  L is  the  first-order  language  of  signature  K,  from 
the  Tarski’s  model-theoretic  truth  definition,  when  a sentence  of  L is  true  in^,  an 
assignment  of  elements  of  A to  variables  satisfies  a formula  of  Z,  in  ^4.  Or  when  the  set  of 
u-tuples  of  elements  of^  that  is  defined  by  a formula  an  «-tuple  (ai,...,a„)  is 

in  the  defined  set  if  and  only  if  the  assignment  taking  each  v,  to  a,  satisfies  the  formula. 

Elementarily  equivalent.  Since  the  signature  K can  be  recovered  from  the  first- 
order  language  L that  it  generates,  the  structures  of  signature  K will  be  referred  to  as  L- 
structures.  Two  L-structures  that  are  models  of  exactly  the  same  sentences  of  L are  said 
to  be  elementarily  equivalent.  Elementary  equivalence  is  an  equivalence  relation  on  the 
class  of  all  L-structures. 

Completeness.  The  set  of  all  the  sentences  of  L that  are  true  in  the  L-structure  A 
is  called  the  complete  theory  of  A,  in  symbols  Th(Y4).  A theory  that  is  Th(v4)  for  some 
structure  A is  said  to  be  complete.  (By  the  completeness  theorem  for  first-order  logic,  a 
theory  is  complete  if  and  only  if  it  is  maximal  syntactically  consistent.)  The  two 
structures  A and  B are  elementarily  equivalent  if  and  only  if  Th(^)  = Th(5). 

Basic  mapping.  The  basic  maps  between  structures  of  the  same  signature  K are 
called  homomorphisms,  defined  as  follows.  A homomorphism  from  structure  A to 
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structure  5 is  a function /from  dom(^)  to  dom(5)  with  the  property  that  for  every  atomic 
formula  ^^(vi,...,v„)  and  any  «-tuple  a = (ai,...,a„)  of  elements  of^, 

^ ^ [a]  =>  5 1=  ^ [b]  L-3 

where  b is  If  we  have  in  place  of  ‘ =>  ’ in  the  quoted  condition,  we  say 

that /is  an  embedding  of  A into  B.  Since  the  language  includes  =,  an  embedding  of  A into 
B is  always  one-to-one,  though  it  need  not  be  onto  the  domain  of  B.  If  it  is  onto,  then  the 
inverse  map  from  dom(5)  to  dom(vl)  is  also  a homomorphism,  and  both  the  embedding 
and  its  inverse  are  said  to  be  isomorphisms.  It  may  be  said  that  two  structures  are 
isomorphic  if  there  is  an  isomorphism  from  one  to  the  other.  Isomorphism  is  an 
equivalence  relation  on  the  class  of  all  structures  of  a fixed  signature  K.  If  two  structures 
are  isomorphic  then  they  share  all  model-theoretic  properties;  in  particular  they  are 
elementarily  equivalent. 

Problem  Statement 

Gallium  nitride  (GaN)  metal  organic  chemical  vapor  deposition  (MOCVD),  from 
the  chemical  engineering  perspective,  is  a convolution  of  several  processes  such  as  heat 
transfer,  species  mass  transfer,  homogeneous  chemical  reactions,  and  heterogeneous 
chemical  reactions  [7].  The  process  may  be  expressed  as 

R^=G{w,T,P,x,}.  1-1 

where  Rg  is  the  growth  rate  of  the  thin  film,  G may  be  view  as  the  operator  representing 
the  MOCVD  process,  y is  the  flow  velocity,  X[  is  the  mole  fraction  of  each  species  and  T 
and  P are  the  system  temperature  and  pressure  . G consists  of  momentum,  energy  and 
species  conservation  operations  and  the  chemical  reaction  operations  and  parameters  1g 
{ju,  Q),  k,  Cp,  Vi,  •••),  where  //  is  the  viscosity  of  the  fluid,  S>is  the  diffusivity,  k is  the 


7 


thermal  conductivity,  Cp  is  the  heat  capacity  and  v,  is  the  stoimatric  parameter  of  each 
species.  G may  be  divided  into  two  categories,  written  in  the  algebra  form, 

G = [F©i?]  1-2 

where  F represents  the  fluid  dynamics,  R represents  the  chemical  reactions  and  © is  the 
convolution  operation.  From  the  transport  phenomena  [5]  and  chemical  reaction  kinetics 
[6],  F and  R are  functions  of  temperature  and  species  concentrations,  and  therefore  are 
self-consistent. 

In  model  theory,  G may  be  seen  as  the  theory  of  the  MOCVD  process  and  1g  iji, 

(D,  k,  Cp,  Vi,  •••)  is  its  interpretation.  Without  1g  (jU,  &,  k,  Cp,  V/,  •••),  G is  only  a set  of 
algebraic  equations  that  states  neither  truth  nor  false.  When  parameters  Ig  (//,  k Cp,  v,, 
etc.)  are  added,  G will  provide  a prediction  of  the  MOCVD  process.  According  to  the 
definition  of  Tarski’s  model  theoretical  truth,  G is  true  if  Ig  {fi,  % k,  Cp,  v„  ■•■)  satisfies 
G.  To  find  Ig  {jn,  % k,  Cp,  v„  •••)  that  satisfies  G is  to  interpret  G.  The  science  of 
chemical  engineering  is  to  perform  the  predictions  and  interpretations. 

Traditional  Chemical  Engineering  Approach 
In  general,  a self-consistent  operation  can  be  deconvoluted  by  trial  and  error,  if 
the  system  is  linear.  When  the  operation  is  nonlinear,  however  multiplets  will  certainly 
destroy  the  uniqueness.  Therefore,  because  of  the  self-consistent  and  high  nonlinearity  of 
operators  F = /(y,  T,  P,  jc, ) and  R = r{T,  P,  x, , M) , deconvolution  of  G will  have  to  have 

the  prior  knowledge  of  three  dimensional  spatially  resolved  gas  phase  temperature  (7), 
pressure  (P),  species  concentration  (xj)  and  reaction  mechanism  (M).  Lacking  such 
information,  the  chemical  engineer  traditionally  would  physically  separate  them.  This 
physical  separation  ensured  the  mapping  of  a chemical  process  into  the  fluid  dynamic 
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“space”  and  the  chemical  reaction  engineering  “space”  where  related  individual  physical 
and  chemical  properties  of  the  process  may  be  studied  separately.  For  example,  under 
certain  process  conditions  such  as  isothermal  ideal  flow,  a deconvolution  can  be 
completed,  by  which  the  original  process  G may  be  mapped  into  the  “space”  chemical 
reaction  engineering,  where  an  individual  process  such  as  homogeneous  ehemical 
reaction  can  be  studied. 

|G  =>  F'  e Fluid  Dynamics  ^ ^ 

[G  =>  i?'  € Reaction  Kinetics 

The  results  of  these  separated  studies  may  then  be  integrated  into  a convoluted 
form,  the  inverse  mapping,  namely  chemical  vapor  deposition  (CVD)  numerical 
simulation  [7-1 1],  in  which  the  entire  CVD  process  including  the  homogeneous  and 
heterogeneous  chemical  reactions,  momentum,  thermal  and  species  conservation 
equations  with  their  boundary  conditions  may  be  solved  simultaneously.  Study  in  this 
separated  form  requires  the  collaboration  of  many  different  research  groups  from 
different  discipline  areas.  It  can  also  be  very  time  consuming  and  often  not  possible  to 
complete.  Worst  of  all,  many  highly  chemical  active  species  just  cannot  be  studied 
separately  (e.g.  reactive  intermediates),  which  will  leave  a high  degree  of  uncertainty  in 
the  final  result. 

\F'^F 

1-4 

[R'^R 

Therefore  often  times,  the  inverse  mapping  does  not  exist,  which  means  that  the  most 
important  tasks  of  chemical  engineering,  the  reactor  design  and  the  process  optimization. 


can  not  be  performed  in  a logically  closed  form. 

g = [f©/?]^  FI®  hi- 


1-5 
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To  overcome  these  shortcomings,  direct  decovolution  is  proposed. 

Approach  of  this  Dissertation 

Mathematical  direct  deconvolution  of  G will  have  to  have  knowledge  of  three- 
dimensional,  spatially  resolved  gas  phase  temperature  (7),  pressure  (P),  species 
concentration  (jCj)  and  reaction  mechanisms  (AT).  Thus,  direct  decovolution  will  have  to 
start  with  a reliable  in  situ  method  that  can  provide  such  information. 

Classical  Theory  of  Raman  Scattering 

Of  all  the  analytical  techniques  applied  to  fluid  inclusions,  only  Raman 
spectroscopy  and  microthermometry  can  provide  information  on  the  composition  and 
pressure/density  of  an  individual  inclusion.  However,  in  microthermomertry,  the 
inference  of  those  two  parameters  is  not  uniquely  determined  from  actual  measurements. 
Therefore,  Raman  spectroscopy  is  the  only  method  that  can  provide  the  information 
needed.  Raman  spectroscopy  has  been  used  in  in  situ  combustion  temperature  studies  for 
decades  [12-14]  and  it  remains  the  only  effective  non-perturbing  way  to  measure  the  gas 
phase  temperature. 

The  phenomenon  of  inelastic  light  scattering  from  a molecule  is  known  as  the 
Raman  effect,  after  the  Indian  physicist  Chandrasekhara  Raman  (1888-1970),  who 
received  a Nobel  Prize  for  the  discovery  in  1930.  Actually,  the  phenomenon  itself  was 
probably  first  recorded  by  R.  W.  Wood,  who  thought  that  the  Raman  lines  were  just 
smudges  on  his  spectrographic  plates  [15].  When  monochromatic  optical  radiation  is 
incident  upon  optically  transparent,  dust-free  gases,  liquids  or  solids,  scattered  light  not 
only  at  the  incident  frequency  (the  process  of  Rayleigh  scattering),  but  also  a very  small 
fraction  of  the  incident  light  scattered  with  a change  in  frequency.  This  was  because  the 
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illuminated  sample  can  absorb  one  photon  while  simultaneously  emitting  another  photon 
of  a different  frequency.  For  a given  molecule  the  induced  dipole  moment 

/u,  = a-E  1-6 

where  a is  the  polarizability  tensor  and  E is  the  external  radiation  field.  Let  E takes  the 


wave  form, 

E = Eq  sm27TUit 


1-7 


then 

IX I = aE  = oEq  sin  Inv,  t 1-8 

Chandrasekhara  Raman  predicted  that  the  polarizability  tensor  may  be  changed  by 
periodic  internal  motion  (rotation  or  vibration  describing  as  an  amplitude  A and 
frequency  v|), 

a =a  + Asm2nVft  1-9 

therefore  the  induced  dipole  moment  becomes, 

Hj  =a  E = {a  + Asm.2nv^t)EQsm27iVit  1-10 

which  may  be  rearranged  as 


Hi  = sin  2kVi  ^ ~ 2k{v,  -i-  u,  )t  +]^AEq  cos  2k{Vi  - u,  )t 


1-11 


These  three  terms  in  turn  are  known  as  the  Rayleigh,  anti-Stokes  and  Stokes  lines  of  a 
Raman  spectrum.  Anti-Stoke  and  the  Stoke  lines  are  called  Raman  shift  and  the  scattered 
light  has  a slight  frequency  shift  due  to  molecular  rotation  or  vibration,  which  provides 
rich  information  of  the  molecules  under  probing.  The  intensity  of  these  shifts  can  be 
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interpreted  through  the  laws  of  physics  to  obtain  rotational  temperature  and  vibrational 
temperature  and  the  molecular  density  of  each  species.  Since  it  is  based  on  the  scattering 
of  light  due  to  an  induced  dipole  moment,  it  can  be  used  to  probe  vibrations  without 
dipole  moments.  Hence,  it  is  applicable  to  H2,  N2  and  O2,  which  are  infrared 
spectroscopic  (IR)  silent. 

Because  of  the  light  scattering  nature  of  spontaneous  Raman  phenomena,  in  a 
typical  Raman  experiment,  a high  intensity  light  beam  is  focused  into  a very  small 
volume.  Therefore,  gas  phase  in  situ  Raman  experiments  can  be  used  to  obtain  three- 
dimensional,  spatially  resolved  gas  phase  temperature  and  species  concentration 
information  to  within  the  limit  of  the  detection  volum.  Gas  phase  temperature 
measurement  using  N2  rotational  transitions  and  species  concentration  measurement 
using  relative  vibrational  transition  intensities  at  room  temperature  are  well  established 
[13,14,16,17].  For  species  concentration  measurement  at  high  temperature,  some 
modification  will  be  needed  and  is  presented  in  Chapter  3. 

Inverse  Problem 

Determine  /f,r(/^,  k Cp,  v;,  •••)  of  F and  R from  three-dimensional,  spatially 
resolved  gas  phase  temperature  and  species  concentration  information  is  an  ‘inverse 
problem’.  It  is  essential  to  find  and  prove  the  uniqueness  in  solving  an  inverse  problem. 
In  general,  the  laws  of  physics  provide  a way  for  computing  the  values  of  the  data  for  a 
given  model.  This  is  called  the  ‘forward  problem’.  In  the  ideal  case  there  is  an  exact 
theory  that  prescribes  how  the  data  should  be  treated  to  reproduce  the  model.  But  in 
realistic  experiments  that  reside  in  a discrete  space,  one  has  a finite  amount  of  data  to 
reconstruct  a model  that  resides  in  a continueous  space  with  infinitely  many  degrees  of 
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freedom  necessarily.  As  the  result,  in  general  the  inverse  problem  is  an  under  determined 
problem,  there  might  be  numerous  models  that  can  explain  the  data  equally  well  yet  can 
not  reproduce  the  truth.  This  is  why  inverse  problems  are  often  called  ill  posed  problems 
[18].  The  illness  here  includes  the  following: 

• Finite  number  of  observations  versus  the  infinite  degree  of  freedom 

• Inherently  noisy  experimental  observations 

The  first  is  related  to  the  uniqueness  of  the  problem  and  the  second  is  related  to  the 
stability  of  the  problem,  which  will  be  discussed  in  detail  in  Chapter  5. 

When  7f,r(/^,  k,  Cp,  v;,  •••)  of  operators  F and  R can  be  uniquely  determined,  the 
mappings  of  convolution  and  the  deconvolution  become  isomorphic  [19]: 

G « [F(y,  T,P,x,)]®[R{T,P,x,,M)]  1-12 

CVD  Modeling 

According  to  Tarski’s  truth  definition,  the  truth  can  not  be  provided  by  G itself 
That  is  a CVD  model  can  not  provide  truth  for  itself 

According  to  the  definition  of  Tarski’s  model  theoretical  truth,  when  the  set  of  n- 
tuples  of  elements  of  A that  is  defined  by  a formula  an  w-tuple  is  in 

the  defined  set  if  and  only  if  the  assignment  taking  each  and  every  v,  to  a,  satisfies  the 
formula.  When  the  conditions  of  the  uniqueness  being  established  and  satisfied,  the 
identification  of  /f,r(a»  2),  k,  Cp,  v„  ••■)  of  operators  F and  R becomes  an  over  determined 
problem  therefore  can  be  solved  uniquely,  Tarski’s  model  theoretical  truth  condition  can 


be  satisfied. 
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When  /f,r(/^,  k,  Cp,  v„  •••)  of  operators  F and  R can  be  uniquely  determined,  G 
and  [i^]©  [r]  are  isomorphic  and  therefore  elementarily  equivalent.  According  to  model 
theory  they  should  share  all  model-theoretic  properties.  Therefore, 

Ig(jU,  % k Cp,  Vi,  -)=/f,r(/^,  S),  k,  Cp,  Vi,-)  1-13 

Therefore, 

For  all  G,  true  if  and  only  if  G’  L-4 

where  G’  is  the  metalanguage  of  G,  which  includes  Equation  1-12  and  7g(//,  k,  Cp,  v,, 
k,  Cp,  Vi,  ■••),  the  unique  solution  of  the  inverse  problem. 

Equation  1-12  will  allow  the  study  of  a CVD  process,  its  reactor  design  and 
process  optimization  in  an  integrated  and  logically  closed  form. 

Overview 

The  overall  plan  of  this  dissertation  is  to  start  with  the  apparatus  and  experiment 
to  examine  MOCVD  of  GaN  by  in  situ  Raman  spectroscopy  (Chapter  2).  Then  the 
measurement  of  gas  phase  temperature  of  species  concentration  is  addressed  and  the 
method  of  data  reduction  discussed  (Chapter  3).  The  next  three  chapters  reported  results 
of  experiments,  TMGa-NH3-N2  gas  phase  reaction  mechanisms  (Chapter  4),  species 
binary  diffusivity  and  gas  phase  reaction  rate  constant  identification  (Chapter  5),  the 
calculation  of  multi-component  system  species  diffusion  and  thermal  diffusion 
coefficient  is  discussed  in  Chapter  6,  while  analysis  of  the  GaN  MOCVD  process 
operating  window  is  presented  in  Chapter  7.  The  information  is  combined  to  present  a 
numerical  simulationof  the  TMGa-NHs-Ni  flow-reaction  system  (Chapter  8),  and  finally 
the  conclusions  of  the  chapters  are  presented  along  with  suggestion  and  future  work 


(Chapter  9). 


CHAPTER  2 

GAS  PHASE  in  situ  RAMAN  SPECTROSCOPY  APPARATUS 
Spontaneous  Raman  Scattering 

In  Chapter  1,  classical  theory  of  Spontaneous  Raman  Scattering  (SRS)  was  briefly 
discussed.  From  a more  quantum  mechanical  point  of  view,  one  may  consider  the  SRS 
effect  to  involve  virtual  states.  In  quantum  mechanical  theory,  the  energy  of  the 
excitation  light  {h  vq)  will  be  absorbed  by  the  probed  molecules,  which  will  be  excited 
onto  a virtual  state.  The  excited  molecules  will  immediately  drop  back  to  their  own  real 
rotational  and  vibrational  states.  The  energy  diagrams  are  shown  in  Figure  2-1 . 

The  vast  majority  of  the  molecules  will  drop  back  to  the  original  rotation  and  vibrational 
states  and  produce  the  Rayleigh  line.  The  molecules  that  drop  back  to  the  same 
vibrational  state  but  different  rotational  states  will  produce  Stokes  (right  wing)  and  anti- 
Stokes  (left  wing)  of  pure  rotational  lines  as  shown  in  Figure  2-2  for  the  spectrum  taken 
for  N2  in  the  experimental  Raman  system  used  in  this  chapter.  The  molecules  that  drop 
to  same  rotational  states  but  one  quantum  number  higher  or  lower  than  the  original 
vibrational  states  will  produce  Stokes  and  anti-Stokes  Q branches  of  the  vibrational  lines, 
as  shown  in  Figure  2-3.  Therefore  the  vibrational  Q branch  is  a collection  of  many 
Raman  transitions.  As  discussed  in  Chapter  3,  the  hot  bands  of  the  SRS  Q branch  are 
produced  by  molecules  that  originated  from  higher  vibrational  states,  will  not  be 
superimposed  onto  the  Q branch  fundamental  band,  molecules  originated  from  ground 
state,  but  produce  a series  of  Q branch  hot  bands.  Vibrational  branches  such  as  P and  R 
branches  are  not  discussed  here. 
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Figure  2-1 . Energy  diagrams  of  of  Spontaneous  Raman  Scattering. 

Instrumentation 

The  essential  elements  of  the  apparatus  necessary  to  carry  out  spontaneous  Raman 
spectroscopy  consists  only  of  a fixed-frequency  laser  source,  a spectrometer  (e.g.  a 
monochromator  to  disperse  the  scattered  light),  and  a detector. 
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Raman  Counts 


Figure  2-2.  N2  Stokes  (right  wing)  and  anti-Stokes  (left  wing)  pure  rotational  lines  taken 
at  room  temperature  with  a Ar^  laser  1.5  W at  488  nm.  Rayleigh  line  was 
filtered  out. 


Raman  Counts 


Figure  2-3.  N2  vibrational  Stokes  line  taked  at  room  temperature  with  a Ar^  laser  6.0  W 
at  488  nm. 
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Light  Source 

SRS  effect  has  one  serious  deficiency:  the  extremely  low  efficiency  of  the 
inelastic  light  scattering.  The  cross  section  for  Raman  scattering  is  typically  lO"'^  of  an 
infrared  absorption  cross  section,  and  the  intensity  of  the  scattered  light  is  1 O'  to  1 0'  of 
the  intensity  of  the  incident  light  [20],  The  low  photon  fluxes  available  from 
conventional  light  sources  meant  that  getting  good  Raman  spectra  requires  exposure 
times  in  the  range  of  minutes  to  many  hours.  Consequently,  spectroscopy  for  real-time 
monitoring  and  minor  species  detection  was  essentially  impossible.  The  radiant  power  of 
Raman  scattering  0ris  given  as  [20], 

0R  Qc  a{Vex)OexE0f^i^~^‘  2-1 

where  a(u«)  is  the  Raman  cross  section  in  [cm  ],  Uex  is  the  excitation  light  line 
frequency,  Eq  is  the  excitation  light  irradiance,  is  the  number  density  in  state  /,  and  the 
exponential  term  is  the  Bolzmann  factor  for  state  /.  This  equation  clearly  shows  that 
Raman  intensity  is  directly  proportional  to  the  source  irradiance  Eq  and  fourth  power  of 
its  frequency. 

The  high  luminosity  of  the  laser  output,  its  narrow  line  width,  and  the  small  size 
and  divergence  of  the  beam  were  all  advantageous  for  Raman  scattering.  While  the 
inefficiency  of  the  inelastic  light  scattering  remained  an  unavoidable  constraint  imposed 
by  nature,  the  laser  reduced  the  severity  of  this  constraint  on  the  application  of  Raman 
spectroscopy.  Since  the  excitation  frequency,  Uex,  can  be  any  frequency,  any  laser  can  be 
used  as  the  excitation  source.  However,  most  Raman  experiments  involve  visible  laser 
sources,  since  they  are  readily  available.  Use  of  ultra  violet  (UV)  laser  requires  special 
care  because  it  might  cause  photochemieal  degradation.  Therefore  the  Spectra-Physics 
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high  power  Argon  ion  (Ar^)  laser  was  selected.  It  delivers  up  to  6.0  W at  487.986  nm. 
This  high  power  density  is  crucial  for  probing  the  gas  phase  temperature  as  high  as  1050 
°C. 

Figure  2-4  shows  the  N2  vibrational  fundamental  Q branch  and  P branch  (similar 
to  the  pure  rotational  peaks)  taken  at  room  temperature.  Among  the  P branch  there  is  a 
small  hot  band  Q branch  peak.  The  very  small  intensity  of  this  peak  when  combined  with 
rotational  spectra  in  Figure  2-4  suggests  that  the  high  power  density  of  the  excitation 
laser  did  not  lead  to  any  laser  heating  in  both  the  rotational  and  vibrational  senses. 

counts  Q branch 


Figure  2-4.  N2  Q branch  hot  band  taken  at  room  temperature  with  a Ar"^  laser  6.0  W at 
488  nm. 

Spectrometer 

Due  to  the  low  intensity  of  the  Raman  scattered  light  compared  to  the  intensity  of 
the  excitation  source,  a complicating  factor  when  selecting  a spectrometer  is  the  variable 
background  signal  in  a Raman  spectrum.  The  farther  the  signals  from  the  Rayleigh-line, 
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the  lower  the  source  background  and  the  better  the  performance  of  a single 
monochromator  system  will  be.  For  visible  excitation  the  Raman  scattered  light 
generally  is  shifted  only  20  to  1 000  A from  the  excitation  wavelength,  so  high  stray  light 
refection  in  the  monochromator  is  essential  for  a successful  measurement.  It  is  necessary 
to  use  a double  or  triple  monochromator  to  totally  disperse  the  scattered  light  in  a Raman 
experiment. 

For  a single  monochromator,  the  resolution  is  determined  by  the  width  of  the 
entrance  aperture,  the  angular  dispersion  of  the  grating  and  the  focal  length  of  the 
spectrometer  optics.  Therefore,  the  slit-width-limited  resolution,  AXs,  neglecting 
diffraction  effects,  may  be  written  as 

AXs^2WI{Da-f)  2-2 

where  W is  the  width  of  the  monochromatic  image  of  the  entrance  slit,/is  the  focal 
length  of  the  spectrometer,  and  Da  is  the  angular  dispersion  of  the  grating.  For  Jobin- 
Yvon  RAMANOR  U-1000  1-m,  double  monochromator  the  resolution  and  slit  width 
calculations  are  more  complicated  and  results  can  be  found  in  reference  [21]. 

When  a CCD  detector  whit  a single  monochromator  is  used,  the  wavelength  range 
that  is  recorded  is  limited  by  the  length  of  the  CCD  and  the  linear  dispersion  of  the 
spectrometer.  For  a given  wavelength  range  recorded  onto  the  CCD,  the  spectral  width 
of  a pixel  is  simply  determined  by  the  dispersion.  For  the  spectral  resolution  of  a CCD- 
spectrometer  system,  if  the  recorded  wavelength  range  is  denoted  by  AX,  and  the  CCD 
has  N pixels  in  the  direction  of  dispersion,  then  the  limiting  spectral  resolution,  SXmin  is 
given  by: 


SA„i„  = 2AA/N 


2-3 
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where  the  factor  of  2 is  included  to  satisfy  the  Nyquist  condition  [22].  D.  A.  Sadler  et  al. 
[23]  pointed  out  that  study  on  photographic  plate  detection  can  be  moved  to  CCD 
detection  by  replacing  the  grain  size  of  the  photographic  emulsion  with  the  pixel 
dimension,  Px. 

Gerard  and  Jacquinot  [22]  suggested  that  for  photographic  plate  detection,  the 
optimum  slit  width  should  be  chosen  so  that  W is  approximately  equal  to  the  grain  size  of 
the  photographic  emulsion.  If  the  entrance  slit  is  narrower,  the  resolution  is  not 
improved,  being  limited  by  SAmm-  Therefore  for  Jobin-Yvon  RAMANOR  U-1000 
equipped  with  a Spectrum  One  CCD  by  Jobin-Yvon  with  AA  chosen  from  50  to  250  cm  ', 
the  actual  number  of  pixels  used  is  740,  pixel  diameter  is  27  pm.  Argon  ion  laser  488  nm 
single  line  as  excitation  source,  which  gives  SAmm=  0.54  cm"'.  On  the  other  hand,  if  the 
slit  is  wider,  the  resolution  is  worsened  without  any  resultant  increase  in  the  illumination 
for  a line  spectrum.  Therefore,  the  optimum  slit  width  at  which  the  spectral  resolution  is 
maximized,  while  simultaneously  satisfying  the  Nyquist  condition,  is  given  by 

W^2  2-4 

Considering  the  aberration  effect,  W should  be  considered  as  the  aberrated  slit  width. 
Therefore,  for  a fixed  wavelength,  the  resolving  power,  Rs=  A / JA  is  maximized  when 
the  aberrated  image  of  the  entrance  aperture  is  matched  to  twice  the  width  of  a CCD 
pixel. 

Detector 

Raman  spectra  are  most  often  obtained  using  a photomultiplier  detector,  with 
scanning  of  the  monochromator.  However,  CCD  detection  has  an  advantage  over  single- 
channel scanning  instruments.  CCD  detectors  can  lower  detection  limits  and  increase 
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sample  throughput  and  quantum  efficiency.  Traditionally,  Raman  systems  have  used 
double  (as  used  in  this  work)  or  triple  monochromators  to  eliminate  source  background, 
and  these  systems  remain  superior  for  source  background  rejection.  However,  a single- 
monochromator with  holographic  filter  systems  has  greater  optical  throughput  than  do 
double  or  triple  monochromators.  A spectrum  containing  256  points  that  takes  T seconds 
to  acquire  with  a single-channel  instrument  can  be  obtained  with  a 256-element  detector 
array  in  T/256s.  In  the  single-channel  instrument,  only  1/256  of  the  light  reaching  the 
exit  of  the  spectrometer  is  available  for  detection,  the  remaining  is  blocked  by  the  exit 
slit.  In  the  CCD  case  all  the  light  reaching  the  focal  plane  of  the  spectrometer  can  be 
detected.  In  practice,  a Raman  spectrum  is  dispersed  along  one  axis  of  the  CCD,  and  the 
spectrometer-slit  dimension  extends  approximately  256  pixels  along  the  other  CCD  axis. 


Signal  Noise  Ratio  of  Photo  Multiply  Tube 

When  PMT  is  used  as  the  detector,  amplifier  readout  and  excess  dark  current 
noise  are  negligible  [20].  Therefore  the  signal  noise  ratio  (S/N)  is  given  as 


where  nR  is  the  Raman  counts,  rif  is  the  number  of  fluorescence  counts  from  the  analyte, 
concomitants  and  the  cell,  rise  is  the  number  of  counts  due  to  elastic  scattering,  rid  is  the 


weak  Raman  signals,  the  S/N  can  be  improved  by 

• increasing  Raman  counts  («/?),  which  can  be  done  by  using  higher  source 
irradiance  or  higher  excitation  frequency, 

• the  number  of  counts  due  to  elastic  scattering  (rise),  which  can  be  done  by  using  a 
double  or  triple  monochromator  that  has  better  stray  radiation  rejection. 


Signal  to  Noise  Ratio  Considerations 


S 


2-5 


number  of  dark  counts,  and  is  the  source  flicker  factor.  Equation  2.5  shows  that  for 
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• the  number  of  dark  counts  {nj),  which  can  be  done  by  cooling  the  PMT, 

• the  source  flicker  factor  (^i),  which  can  be  done  by  using  more  stable  light  source 
such  as  solid  state  laser, 

• and  by  reducing  the  number  of  fluorescence  counts  from  the  analyte, 
concomitants  and  the  cell  {np),  which  is  the  most  difficult  tasks  in  Raman 
spectrometry.  One  solution  is  to  change  the  excitation  frequency.  But  since 
Raman  counts  is  proportional  to  the  fourth  power  of  the  excitation  frequency. 
Reducing  excitation  frequency  may  compromise  the  Raman  counts. 

In  high  temperature  application,  thermal  radiation  may  produce  a high  level  of 

noise.  If  the  temperature  of  the  heater  is  less  than  1 500  °C,  most  thermal  radiation  is  in 

the  inferred  region.  S/N  may  thus  be  improved  by  using  high  excitation  frequency. 

Signal  to  Noise  Ratio  of  Charge  Couple  Device  Array 

In  Raman  signal  the  major  noise  source  allows  the  analyst  to  determine  whether 

binning  or  co-addition  should  be  used. 

Limit  by  Background  Shot  Noise 

Background  photons  are  a major  source  of  noise  in  the  Raman  measurement. 
Based  on  the  limiting  expressions,  the  maximum  S/N  under  these  conditions  is  obtained 
when  background  photon  shot  noise  is  the  dominant  noise  source.  Also  if  a measurement 
is  background  photon  shot  noise  limited,  there  is  no  S/N  difference  between  reading  from 
a single  large  well  and  co-addition  of  data  read  from  multiple  wells.  However,  if  binning 
is  not  used,  cosmic  ray  removal  is  simplified  and  device  full  well  capacity  is  effectively 
increased.  Binning  the  device  will  result  in  a well  capacity  limited  by  the  capacity  of  the 
output  register  and  output  amplifier.  On  the  other  hand,  co-addition  of  the  signal  from 
256  pixels  results  in  a maximum  obtainable  S/N  of  eight  times  what  is  possible  with 
binning. 
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Limit  by  Readout  Noise 

At  liquid  N2  temperature  the  dark  current  of  individual  pixels  should  be  very  low. 
Therefore  for  pure  N2  gas  temperature  measurement,  background  noise  is  also  low. 
Readout  noise  can  be  significant.  A number  of  pixels  corresponding  to  the  slit  image 
height  may  be  binned  together,  resulting  in  a signal  equal  to  the  sum  of  the  signal  in  the 
individual  pixels,  but  with  a readout  noise  equivalent  to  a single  pixel.  The  binning 
process,  by  reducing  the  number  of  readout  operation,  can  significantly  improve  the 
signal-to  noise  ratio. 


Figure  2-5.  Schematic  diagram  of  the  MOCVD  system  with  a in  situ  gas  phase  Raman 
probe. 


24 


A similar  experimental  apparatus  has  been  used  for  previous  studies  [24] . The 
system  has  been  completely  rebuilt  along  with  some  significant  improvements  for 
quantitative  interpretation  purposes.  It  has  a 1995  model  Jobin-Yvon  RAMANOR  U- 
1000,  equiped  with  a double  1-m  monochromator  with  two  plane  holographic  gratings 
(1800  grooves/mm).  The  Raman  scattered  lines  can  be  detected  with  either  a R928-01 15- 
0381  photomultiplier  by  Products  for  Research,  Inc  or  Spectrum  One  CCD  by  Jobin- 
Yvon.  The  Coherent  Innova  90-5  Argon  ion  laser  was  replaced  by  a much  powerful 
Spectra-Physics  high  power  Argon  ion  laser,  which  delivers  up  to  6.0  W at  487.986  nm. 
The  signal  collection  mirror  in  the  macro  chamber  was  taken  out.  It  was  found  that  this 
mirror  would  reflect  the  emission  of  the  heater  to  the  probing  spot,  and  consequently 
reduce  the  signal  noise  ratio.  Since  the  signal  to  noise  ratio,  not  the  signal  intensity 
alone  required  enhancement,  it  was  believed  to  be  better  to  remove  it.  A schematic  of  the 
in  situ  Raman  spectrometer  system  is  shown  in  Figure  2-5.  Spectral  resolution  was 
maintained  at  better  than  2.0  cm'*.  All  calibrations  were  done  by  manufacturer  service. 

Inverted  Impinging  Flow  Reactor 

The  schematic  diagram  of  the  reaction  system  is  shown  in  Figure  2-5.  The  detail 
schematic  diagram  of  the  reactor  is  shown  in  Figure  2-6.  The  reactor  design  provides 
three  separated  circular  gas  inlet  streams.  The  most  inner  center  stream  contains  the  Ga 
source,  trimethylgllium  (TMGa),  while  the  outer  center  stream  contains  NH3.  The  height 
of  the  inner  most  center  tubing  can  be  adjusted  to  achieve  a different  degree  of  mixing. 
The  sweeping  stream  contains  carrier  gas  only  (N2  or  H2),  and  is  designed  to  reduce  the 
concentration  of  reactants  near  reactor  walls  and  therefore  keep  the  optical  windows 
clean  and  at  same  time  reduce  the  re-circulation  generated  at  the  edge  of  the  susceptor  in 


cold  wall  vertical  reactor  [25]. 
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Figure  2-6.  Schematic  diagram  of  a vertical  upflow  reactor. 
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Figure  2-7.  Schematic  diagram  of  reference  position  determaination  setup. 
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Figure  2-8.  Pure  N2  vibrational  Q branch  intensities  along  the  centerline  of  the  reactor. 
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Reference  Point  Position  Determaination 

The  schematic  of  the  reference  position  determaination  procedure  is  demonstrated 
in  Figure  2-7.  Earlier  study  [24]  suggested  that  using  the  position  where  the  N2  intensity 
equals  95%  of  that  from  the  position  at  infinity  as  the  reference  point.  A large 
discrepancy  was  found  when  the  experimental  gas  phase  temperature  results  were 
compared  with  simulated  results.  This  discrepancy  was  thought  [24]  to  be  caused  by 
possible  poor  contact  between  heater  and  the  heater  envelope. 

From  Figure  2-7,  one  may  find  that  at  the  ideal  condition,  the  95%  intensity 
position  should  be  around  5 mm  away  from  the  suscepter  not  0.3  mm  that  the  earlier 
study  assumed.  Because  the  heater  envelope  will  block  the  Raman  signal  solid  angle  first 
not  the  Ar^  laser  beam  solid  angle,  since  Raman  signal  has  much  larger  solid  angle. 

Travling  by  approximately  the  same  distance  from  the  macro  chamber  wall  to  the 
focus  point,  the  Ar"^  laser  beam  reduces  its  diameter  from  8 mm  to  about  0.02  mm,  while 
the  Raman  siginal  reduce  its  irradiant  diameter  from  44mm  (effective  diameter  of  the 
siginal  collection  lens)  to  about  1 mm.  Solid  angle  is  in  cone  shape  as  shown  in  Figure  2- 
7 assuming  the  diameter  at  the  focal  point  is  zero.  When  heater  envelope  moves 
downward,  it  comes  into  contact  with  the  edge  of  the  cone.  Larger  cone  comes  into 
contact  first. 

From  Figure  2-7,  the  radias  of  the  heater  envelope  is  about  20  mm.  The  focal 
length  of  both  Raman  siginal  collection  lens  and  laser  beam  focal  lens  are  the  same  of  80 
mm.  Therefore,  when  the  heater  envelope  moves  downward,  first  in  contact  with  the 
solid  angle  cone  is  the  edge  of  the  heater  envelope,  which  is  20  mm  from  the  focal  point 
at  one  forth  of  the  focal  length.  The  cone  diameter  of  the  Raman  signal  is  about  1 1 mm, 
while  the  cone  diameter  of  the  laser  beam  is  about  2 mm.  Since  Raman  signal  has  a 
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larger  solid  angle  cone  diameter,  it  should  be  blocked  by  the  heater  envelope  first. 
Therefore,  the  first  5%  of  the  signal  loss  certainly  comes  from  the  the  blocking  of  the 
Raman  signal  solid  angle  at  5.5  mm  from  the  susceptor.  If  the  heater  envelope  continues 
to  move  downward,  at  about  1 mm  from  the  the  susceptor,  the  heater  envelope  starts  to 
come  in  contact  with  the  laser  beam  cone.  Figure  2-8  shows  the  experimental  results  of 
the  pure  N2  vibrational  Q branch  intensity  as  a function  of  position  along  the  centerline  of 
the  reactor.  From  this  figure,  one  may  find  that  the  Raman  intensity  starts  to  drop  at 
about  5.5  mm  and  then  drops  dramatically  at  1 mm  from  the  susceptor. 

In  the  earlier  study  smaller  power  Ar^  laser  was  used  [24]  which  has  a beam 
diameter  of  about  3 mm,  one  fourth  of  3 mm  is  0.75  mm.  The  point  that  the  heater 
envelope  starts  to  block  the  solid  angle  of  that  laser  beam  should  be  at  0.375  mm  from  the 
susceptor. 

In  addition  to  the  heater  envelope  blocking  the  Raman  signal  solid  angle,  both  the 
curvature  of  the  cylindrical  quartz  reator  and  its  slit  opening  can  also  influce  the  Raman 
intensity  and  are  yet  difficult  to  quantify.  Therefore,  95%  Raman  intensity  point  should 
be  a point  in  between  5.5  mm  to  1 mm  from  the  susceptor.  However,  there  is  no 
definitive  way  to  determain  exactly  where. 

Therefore,  this  work  chose  the  point  that  the  heater  envelope  starts  to  block  the 
solid  angle  of  the  laser  beam  as  the  reference  point.  The  exact  position  is  calculated  from 
the  beam  diameter  of  the  laser  beam  using  the  method  demonstrated  in  Figure  2-7,  which 
should  give  much  better  reliability  and  repeatbility. 


CHAPTER  3 

GAS  PHASE  IN  SITU  RAMAN  TEMPERATURE 
AND  SPECIES  CONCENTRATION  MESEARMENT 


N?  Rotational  Raman  Spectra  for  Gas  Phase  Temperature  Measurement 
The  energy  of  a classical  rigid  rotator  can  be  found  from  any  classical  mechanics 
text  book,  that  is, 

E = -Ico^  =J^/2I  3-1 

2 

where  J is  the  angular  momentum  and  l=juR^  is  the  rotational  inertia,  fj.  is  the  reduced 
mass  and  R is  the  length  of  the  rotator.  The  quantization  of  angular  momentum  is  given 
by, 

+ 3-2 

where  j is  the  quantum  number,  h is  the  Plank’s  constant.  Then  the  energy  levels  of  the 
quantized  rotator  are 

=[*V(8a-V)l/X;  + l)  3-3 

For  Raman  scattering  the  angular  momentum  selection  rule  follows  two  successive 
electric  dipole  transitions  Aj  = 0,  ± 2.  Aj  = 0,  is  the  Rayleigh  peak,  and  the  pure  rotational 
Raman  Stoke-lines  are  then 

Au  = AEj/hc  = {Ej^^ -Ej)!hc 

= [/2/(8;r'c/)](4y  + 6)  3-4 

= 5(47  + 6) 

where  B is  the  rotational  constant,  signature  of  the  molecule  [28]. 
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The  photon  scattering  rate  Ij  of  any  Raman  line,  which  involves  the  initial 
rotational  state  j,  is  simply  proportional  to  the  population  of  this  state  as  described  by  a 
Boltzmann  distribution.  Thus  the  photon  scattering  rate  is 


I v]  + \){2j  + l)exp 


zh" 

V J 


^ ^ ^as  (Ac,  +l)(2;  + l)exp| 


-E  ^ 
V ^B^  J 


(Stokes) 


, (anti-Stokes) 


3-5 


3-6 


where  Us  and  L>as  are  Stokes  and  anti-Stokes  lines.  Ej  and  Ej+2  should  not  be  confused  with 
the  Raman  shift  he  zL)  The  exponential  term  represents  the  Boltzmaim  distribution  of  the 
number  of  molecules  occupying  on  different  energy  levels  at  a given  temperature.  The 
factor  is  used  for  photon  counting  while  factor  is  used  for  intensities. 

For  a real  molecule  such  as  N2,  however,  the  rigid  rotator  model  is  not  an  accurate 
model,  many  correction  parameters  are  needed  to  modify  the  rigid  rotator  model  [13]. 
After  rearrangement,  the  following  equation  may  be  used  to  calculate  the  N2  rotational 


temperature. 


In 


j<-j+2 


i^o  -svaYs{j)f{j)g. 
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3-7 

where  is  the  incident  laser  frequency  488  nm,  is  the  wave  number  shift  associated 

with  a pure  rotational  transition  {j<—j+2),  and  s=+\  for  a Stoke  shift  and  s =-l  for  a anti- 
Stoke  shift.  The  nuclear  spin  degeneracy  gj=3  for  N2.  The  constant  h and  kQ  are  the 
Planck  and  Boltzmann  constants,  respectively.  Cexp  is  an  experimental  constant  that 
accounts  for  the  efficiency  of  the  detection  system,  and  ^ Pq  ) is  the  spectral 
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response  of  the  spectrometer  and  detection  system.  The  term  f(j)  is  a correction  term  for 
the  anisotropy  of  polarizability  tensor,  (y)o^,  and  it  accounts  for  the  effect  of  molecular 
non-rigidity.  For  N2  this  function  is  given  by 

/0')=l+4z“'(Be/ffle)^b'0'+l)+a+2)0'+3)|  3-8 


where  /=0.45  and  BJco^  = 8.476x1 0'"*  for  N2  [26].  In  Equation  3-7,  S(j)  is  the  rotational 
Raman  line  strength,  which  is: 


S{j)  = 


3(7  + !)(;•  + 2) 

2(27  + 3) 


3-9 


Ej{j)  is  the  energy  of  the  rotational  level  as  give  by 

Ej  U)  = BJU  + 1)  - D,f  (7  + 1)^  + (7  + 1)^  3-10 

The  N2  rotational  constants  were  experimentally  determined  [27].  Bo  =1 .989574  cm'*, 
Z)fl=5. 76x10'^  cm'',  and  Ho  =0.0  cm''. 


= E(j  + 2)  - E(j) 

= (45.  -6D.  «„)(;  + |)-(8D. + + 

3-11 

The  diatomic,  molecular  rotational  partition  function  Qr  is  given  as, 

Qr  = TjSjiV  + ^)QW(-hcF{j)/kgT)  3-12 

J 

From  Equation  3-7  it  is  seen  that  plotting  ln{7,.y3-2/[5'(/’)/(/')gy]}  against  F(j)hclkB,  should 
produce  a straight  line.  The  reciprocal  of  its  slope  equals  T,  the  nitrogen  rotational 
temperature. 

Rotational  temperature  some  times  may  not  be  the  same  as  the  gas  temperature, 
mainly  due  to  the  deviation  from  the  Bolzmann  distribution  caused  by  high  energy 
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sources  (such  as  chemical  or  atomic  reactions,  collision  by  high  energy  particles  etc.). 
Typical  relaxation  time  for  rotational  temperature  is  less  than  10'^  seconds.  The  time  for 
gas  phase  molecules  pass  through  probing  area  is  about  0.5~1  .Ox  1 0'^  seconds.  Therefore, 
the  rotational  temperature  we  measured  should  represent  the  gas  temperature. 

MOCVD  Reactor  Transient  Temperature  Study 

Introduction 

Compound  semiconductor  materials  are  used  in  a wide  range  of  technological 
applications,  including  the  fabrication  of  microelectronic  circuits,  optical  and  magnetic 
recording  media,  optical  devices  and  high-performance  cutting  and  grinding  tools.  For 
this  wide  range  of  applications,  heteroepitaxy  technology  is  often  used,  in  which  a thin 
film  is  grown  on  a substrate  of  different  materials.  For  example  GaN  has  been  grown  on 
sapphire,  silicon,  silicon-carbide  (6H-SiC),  LiGa02  (LGO)  and  LiAlOi  (LAO).  The 
properties  of  the  GaN  grown  on  different  substrate  are  found  to  be  very  different. 
Modeling  results  discussed  in  a later  chapter  indicates  that  the  gas  phase  temperature 
plays  a very  important  role  in  determine  the  gas  phase  reaction  path.  Because  of  the 
differences  in  termal  properties,  the  growth  surface  temperature  dynamics  may  vary  with 
substrate  selection,  Raman  temperature  measurement  provides  an  excellent  method  to 
probe  this  effect. 

Experimental 

An  inverted  stagnant  flow  reactor  with  in  situ  Raman  spectroscopy  was  used  as 
shown  in  Figures  2-5  and  2-6.  A pure  Nitrogen  (99.999%)  gas  flow  was  introduced  at  the 
bottom  of  the  reactor  and  directed  toward  a substrate.  The  gas  velocity  was  set  to  be  at 
3.0  cm/s  based  on  22  °C  and  1 atm  for  all  experiments.  In  this  study,  2’  silicon  and 
sapphire  wafers  were  used  as  substrates  and  one  experiment  was  performed  without  any 
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substrate  as  comparison.  After  the  flow  pattern  reached  a steady  state  at  room 
temperature,  the  substrate  was  heated  by  an  electrical  heater  placed  in  a heater  envelope 
with  a specified  electrical  power.  Temperature  of  the  heater  was  measured  with  a 
thermocouple  placed  in  the  middle  of  the  heater  touching  the  inside  surface  of  the  heater 
envelope.  An  example  recorded  temperature  profile  is  shown  in  Figure  3-1 . Each 
experiment  was  started  at  22  °C,  and  the  maximum  temperature  was  set  at  790  °C. 


Figure  3-1 . Susceptor  temperature  ramp  profile  for  all  transient  experiments. 

In  situ  Raman  spectra  were  taken  at  a fixed  spatial  point  2 mm  from  the  substrate 
in  the  center  of  the  reactor.  The  488  nm  line  of  an  Ar^  laser  at  1 .6  W was  used  as  the 
excitation  source,  and  scattered  Raman  signal  was  collected  at  a right  angle  with  a CCD 
detector  through  a single  monochrometer.  A total  of  30  cycles  were  taken  for  each 
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experiment,  with  an  interval  of  20  seconds  between  each  cycle.  Each  cycle  consists  of 
two  scasn  with  an  integration  time  of  2.5  seconds,  and  the  same  for  the  dark  cycles. 

Counts 
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Figure  3-2.  N2  Raman  spectra  taken  at  5 mm  from  the  susceptor  (no  substrate  was 
present)  by  CCD. 

Results  and  Discussion 

Transient  behavior  of  heat  transfer  in  the  presence  of  different  substrates 
(sapphire  and  silicon)  was  investigated  using  in  situ  Raman  spectroscopy.  N2  rotational 
Raman  spectra  in  the  case  of  no  substrate  presence  are  shown  in  Figure  3-2.  Because 
single  monochrometer  was  used  in  conjunction  with  a CCD  detector,  the  stray  light  was 
not  totally  filtered  out.  The  large  tail  of  the  Rayleigh  line  was  clearly  present  ans  is  seen 
to  interfere  with  the  lower  part  of  the  spectra.  Cosmos  particles  also  cause  large 
interference  in  the  spectra  resulting  loss  of  usable  information  as  shown  in  Figure  3-2. 

The  gas  phase  transient  temperature  profiles  at  2 mm  from  the  substrate  in  the 
center  of  the  reactor  were  calculated  from  in  situ  Raman  N2  rotational  spectra  using 
Eguation  3-7.  The  calculated  results  of  each  spectrum  were  plotted  as  shown  in  Figure  3- 
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3 to  obtain  the  gas  phase  temperature  of  each  spectrum.  The  gas  phase  transient 
temperature  profiles  at  2 mm  from  the  substrate  in  the  center  of  the  reactor  for  cases  of 
bare  susceptor,  sapphire  and  silicon  as  substrate  are  shown  in  Figure  3-4. 


Figure  3-3.  Example  of  the  linear  regression  fit  of  Equation  3-7  for  gas  phase 
temperature  calculation  from  the  N2  rotational  Raman  spectrum. 


The  transient  behavior  of  the  gas  phase  temperature  was  found  to  be  identical  for 
bare  susceptor  and  silicon  cases  at  a high  temperature  ramp  rate.  When  a sapphire  wafer 
was  placed  on,  the  gas  temperature  rise  was  slightly  slower  as  a result  of  its  lower 
thermal  conductivity  and  high  transparency.  In  three  separate  cases,  sapphire  substrates 
broke  into  small  pieces  due  to  thermal  stress  from  the  high  temperature  ramp  rate. 
Therefore  only  uncompleted  experimental  results  were  able  to  be  plotted  in  Figure  3-4  for 


sapphire  case. 
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Figure  3-4.  Gas  phase  transient  temperature  profiles. 

Conclusion 

Despite  interferences  from  different  sources,  the  CCD  detector  produces  a 
reasonable  S/N  for  fast  transient  study  of  the  heat  transfer  behavior  of  the  gas  phase 
during  heater  ramp-up.  The  transient  behavior  of  the  gas  phase  temperature  was  found  to 
be  similar  for  bare  susceptor  and  silicon  cases  at  very  high  temperature  ramp  rate.  For 
sapphire  the  rate  of  the  gas  phase  temperature  rise  was  slightly  slower. 

Vibrational  Raman  Intensity  for  Species  Concentration  Measurement 
Relative  Normalized  Differential  Raman  Scattering  Cross-Section 

The  Intensity  of  a vibrational  Raman  transition  is  directly  proportional  to  the 
number  of  molecules  as  Equation  3-13 


where  Ij  is  the  vibrational  intensity  of  the  scattering  molecule,  /o  is  the  power  density 
of  the  incident  radiation,  Nj  is  the  number  of  scattering  molecules,  and  (5cr  / dQ)  ^ as 

given  in  Equation  3-14  [16]is  the  differential  scattering  cross-section  for  the /*’ 
vibrational  band  and 

=^(S«  -^.)‘(v,  +n/y;g,[45(a)f+7(r,)]  3-14 

5Q  45  j j j 

in  this  equation  vq  is  the  excitation  line  source  frequency  [cm''],  Vj  is  the  Stokes  shift  of 
the  f'  vibrational  peak  frequency  of  the  scattering  molecule  [cm''],  £>/  =hl^7^cVj  is  the 
square  of  the  zero  point  amplitude,  h is  the  Planck's  constant  [28],  c is  the  velocity  of 
light  in  vacuum,  g,  is  the  degree  of  degeneracy  of  the  vibrational  state  [29]  and  a],  and  yj, 
are  the  components  of  polarizability  tensor  [30], 

The  measurement  of  species  abundance  can  be  made  by  measuring  scattered 
intensity  using  the  relation  in  Equation  3-15  by  given  the  absolute  Raman  differential 
cross  sections. 

3^5 

{daldO),^ 

A direct  determination  of  the  absolute  Raman  differential  cross-section, 

(dcr  / 5Q)  p however,  is  often  difficult  and  tedious  [16].  Therefore,  the  analysis  of 

scattering  intensities  is  normally  made  on  a relative  basis.  In  this  approach  the  scattering 
intensity  of  one  species  of  unknown  concentration  in  a mixture  is  measured  only  with  that 
of  a reference  species  in  the  same  mixture  of  known  quantities  and  a well  established 
differencial  scattering  cross-section.  This  concept  has  been  explored  by  Weber  and 
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Brodersen  [16],  who  selected  the  N2  Q-branch  vibrational  transition  at  2331  cm-1  as  the 
reference  signal.  In  their  analysis,  a "Relative  normalized  differential  Raman  scattering 
cross-section  term  Sj  was  introduced  and  defined  bythe  following  equation:  Numerous 
gaseous  species  were  determined  and  tabulated  for  several  excitation  lines. 


2,  = 


(acr/ao), 


(uo  -Oj) 


-4 


' (ao-/aQ)o.  (uo -233 lew-') 


1 ^-4 


• [1  - QX^{-V]hc  / kT)\ 


2331cm-'  gjma'if+lXjirjY} 


This  formular  has  been  the  basis  for  quantitative  measurement  of  gas  phase  composation 
by  Raman  spectrum  as  reviewed  in  the  literature  [17,31].  It , however,  that  the 
found  it  is  quite  necessary  to  do  some  modification  in  order  to  extend  this  concept  to  high 
temperature  gas  phase  species  abundance  measurement. 

Modification  of  the  Relative  Normalized  Differential  Raman  Scattering  Cross  Section 
The  relative  normalized  differential  Raman  scattering  cross  section  formula  of 
Equation  3-16  was  designed  for  concentration  calculation  at  room  temperature.  In  order 
to  calculate  concentration  at  elevated  temperature,  a new  formula  needs  to  be  established. 
Simple  Modification 

The  primary  modification  is  to  account  for  the  change  in  the  population  of 
fundamental  state  of  N2.  Adding  the  temperature  correction  term  to  the  nitrogen 
intensity,  gives  the  following: 


, ^ (5cr /5Q)^  {vq-Vj)^  \-QXTp{-Vjhc I kT) 

{6a/dQ.)Qjj^  (l>o -2331cm-')-'  l-exp(-2331cm-'/ic/A:r) 
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therefore 

E'=E, j 3-18 

^ 1 -exp(-2331cA?j  ^hcIkT) 

Critical  Review  of  Vibrational  Raman  Intensity 

Raman  vibrational  spectra  of  N2  at  elevated  temperature  for  the  purpose  of 
vibrational  temperature  measurement  were  reported  by  Lapp,  Glodman  and  Penney  [12], 
The  spectra  indicate  that  the  assumption  that  “transition  l<-2,  2<— 3,  •••  will  overlap  with 
the  0<-l,  as  the  fundamental  bands”  is  not  entire  correct.  Instead,  hot  bands  Raman 
transitions  red  shifted  as  much  as  26  cm’’  towards  the  excitation  line.  This  is  because  the 
vibrational  energy  of  N2  is  given  as 

G{U,  J)  = E{V,  J)/  hc  = COe{p  + 1 / 2)  - COeXeiV  + 1/2)^  + tvt;ye(o  + 1/2)^ 

+ (Be-Oe/2)J(J  + l)-(De  + J0c/2)j\j  + iy  - C(eU/(J  + 1)  ~ (J  + if  + ' ' ■ 

3-19 

where  cOe,  (o^e,  and  cOeye  are  vibrational  constants;  Be  and  De  are  the  rotational  constants; 
and  Oe  and  are  the  vibration-rotation  interactions.  The  Raman  shift  according  to  the 
selection  rule  for  the  Q-branch  is  Av=  1 and  AJ=0.  Hence,  the  N2  hot  band  red  shift 
will  be 

AG(u  + \,J  <^u,J)-cOe-2o}eXe{u  + \)  + cOeye(3u^  +6lh-13/4)  ^ 

It  is  clear  that  this  red  shift  originates  from  the  anharmonic  behavior  of  the 
vibrational  energy  potential.  Since  AG  is  a function  of  v,  they  can  not  overlap  onto  each 
other.  For  the  Raman  spectrometer  system  usually  show  a spectral  resolution  no  worse 
than  4 ~5  cm’’.  The  rransitions  from  hot  bands  should  be  clearly  resolved  since  they  are 
separated  at  26  cm’’.  The  N2  Q-branch  2331  cm’’  represents  only  the  transition  from 
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ground  population  at  both  high  and  low  temperature. 


The  same  argument  should  apply  to  the  vibrational  Raman  intensity  of  other 
species.  For  example,  ammonia  has  four  Raman  active  vibrational  Raman  transitions. 


Table  3-1  listed  their  v=l  and  v=2  vibrational  level  frequencies  and  related  Stokes 


Raman  shifts.  These  data  clearly  show  their  anharmonicity  and  therefore  should  not 


overlap  with  each  other. 


Table  3-1.  Ammonia  vibrational  frequencies  and  related  Stokes  Raman  shifts  [32] 


vl 

v2 

v3 

v4 

V=1 

3336.0  3337.1 

932.4  968.1 

3443.6  3443.9 

1626.1 

1627.4 

v=2 

1598.5  1882.2 

6850.0  6850.4 

3216.1 

3217.8 

Raman  0<—  1 

3336.0  3337.1 

932.4  968.1 

3443.6  3443.9 

1626.1 

1627.4 

(Stokes)  l<-2 

3278.3  3278.9* 

666.1  914.1 

3406.4  3406.5 

1590.0 

1590.4 

* No  experimental  data  available  for  high  level  vibrational  states  of  vl  mode,  instead, 
theoretical  calculation  results  are  used  here  [33]. 


Intensity  of  Vibrational  Raman  Transition  from  0<— 1 Only 

Since  Raman  intensity  is  proportional  to  the  population  Nj  that  participates  in  the 


given  Raman  transition  (Equation  3-13),  we  need  to  use  the  relationship  between  only  the 


ground  population  and  the  total  population  of  the  species  N in  the  sample  for  the 


cross-section  should  be  used.  Retaining  the  Boltzmann  expression  for  the  population 
in  different  states 

J 


A'  = A^v,=oO-exp‘-'""'^>)-' 


3-21 


and  Vj=0, 


o'*  .„-'* 

L n — 


45 


(uo  -vj)\\-e 


(-hif/kT) 


)b^g\45(a)r+lzXry] 
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The  ratio  of  the  number  of  0<— 1 transition  out  of  a total  v<— v+1  transitions  is 
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3-23 

With  this  expression,  the  modified  "Relative  normalized  differential  Raman  scattering 
cross  section"  is  as  follow: 

^ {d(jldO)j  {uo-Ojy'^  \ -exp(-233 \cm~^ he / kT) 

^ (Scr/SQ)^  ^^  (uo -233 1cm 1 -exp(-u,7ic/A:r) 

_ 2331cm~'  +lXjirjY}  ^ 24 

mak.Y+ix.Srkf] 

The  conversion  factor  is 

_ 1 - exp(-2331cm’'/ic/A:T)  2 25 

[l  - exp(-3337cm-'/?c/itr)f 


Table  3-2.  NH3  and  TMGa  mole  fractions  from  2~,  27  and  I” 


Temperature(K)  NH3  mole  fractions*  TMGa  mole  fractions** 


s 

I’ 

S” 

E 

E’ 

E” 

300 

0.227 

0.227 

0.227 

0.002282 

0.002282 

0.001932 

500 

0.227 

0.227 

0.227 

0.002689 

0.002686 

0.001641 

700 

0.227 

0.226 

0.228 

0.003172 

0.003146 

0.001401 

900 

0.228 

0.224 

0.231 

0.003683 

0.003595 

0.001225 

1100 

0.229 

0.221 

0.234 

1300 

0.232 

0.218 

0.237 

1500 

0.235 

0.215 

0.240 

* Calculated  from  relative  Raman  intensity  of 


iNi 


= 0.5 


** 


Calculated  from  relative  Raman  intensity  of 


I TMGa 

In, 


0.001 


The  vibrational  Raman  cross-section  (5cr  / 5Q)  , of  TMGa  reported  by 

526cm~\TMGa  ^ 


[34]  is  1.3x10  ^^ern^sr  ' per-molecule  at  28  °C  for  488  nm  excitation.  The  reported 
vibrational  Raman  cross-section  (5cr/5Q)g^^  ofN2  [16]  is: 

5.05  X - 2331cm~  ^ at  28  °C  for  488  nm  excitation.  Therefore 
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the  modified  relative  normalized  differential  Raman  scattering  cross  section  of  TMGa 

ti 

S is  17.496  at  room  temperature.  The  detail  of  the  relative  eoneentration 

526cm  , TMGa 

method  and  the  modified  Relative  Normalized  Differential  Raman  Scattering  Cross 
Section  formula  can  be  found  in  [35]  and  the  references  therein. 

The  caleulated  NH3  and  TMGa  mole  fractions  from  1]  F and  F”  are  listed  in 
Table  3-2  for  comparison. 


CHAPTER  4 

GaN  MOCVD  GAS  PHASE  REACTION  MECHANISM  STUDY 


Device  quality  Gallium  Nitride  (GaN)  thin  films  sere  finally  demonstrated  in  the 
early  1990s,  largely  due  to  a innovative  MOCVD  reactor  design  [36,37],  Since  then 
several  studies  have  been  dedicated  towards  understanding  the  reaction  mechanism  [38- 
41],  After  more  than  one  decade,  the  mechanism  is  still  not  fully  understood. 

Introduction 

TMGa  Thermal  Decomposition 

The  thermal  decomposition  of  TMGa  is  the  most  extensively  studied  species,  as 
both  ab  initio  quantum  mechanical  calculation  and  experimental  studies  have  been 
repeated  by  several  authors.  A gas  phase  decomposition  mechanism  was  also  proposed 
[42-45]  and  reaction  rate  constants  for  each  step  have  been  reported  [42,43].  It  is  clear 
that  trimethylgallium  (TMGa)  decompose  to  form  dimethylgallium  (DMGa)  then 
monomethylgallium  (MMGa), 


ki 

TMGa  DMGa 


k2  ks 

^ MMGa  ^ Ga  R-1 


The  reported  rate  constants  are  summarized  in  Table  4-1. 


Table  4-1 . TMGa  pyrolysis  rate  constants 

ki  ko  [s‘‘]  3.5xl0‘^ 

Ea  [kcal]  59.5  64.89 

k2  ko[s'‘]  8.7x10^ 

Ea  [kcal]  35.4 

Ks  ko[s'']  l.Oxlo'^ 

Ea[kcal]  77.5 

References  [42]  (1962)  [44]  (1988) 


2.9xl0‘^ 

64.6 

4.0x10'^ 

52.6 
l.OxlO'^ 
54.1 

[43] (1991) 


2.37x10'^ 

47.1 


[45]  (1994) 
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The  TMGa-N2  system  reaction  rate  constant  ki  was  identified  and  will  be 
discussed  in  Chapter  5. 

NHi  Homogeneous  Decomposition 

The  NH3-Ar-system  has  been  studied  by  several  groups  using  the  so-called  shock 
tube  method.  These  results  reported  a slow  second  order  reaction  for  the  initial  step 
followed  by  many  secondary  reactions: 

ki 

NH3  + M ^ NH2  + H R-2 

Values  of  the  rate  constant  of  this  initial  step  from  different  groups  are  summarized  in 
Table  4-2. 

Table  4.2  NH3  pyrolysis  rate  constant  [cm^  mof‘  s'*] 

ki  ko  [cm^  mof'  s'']  1.2x10*^  4.0x10*^  l.gxio'^  2.2xl0'" 

Ea[kcal]  91.0  93.9  92.1  93.5 

Reference  [46]  (1979)  [47]  (1981)  [48]  (1981)  [49]  (1990) 


The  NH3-N2  system  reaction  rate  constant  ki  was  identified  and  will  be  discussed 
in  Chapter  5. 

TMGa-NH^ 

Several  research  groups  have  studied  the  gas  phase  the  reactions  of  TMGa-NH3- 
H2  system  [38,39  ,50  -53].  It  is  widely  accepted  that  at  low  temperature,  especially  at 
room  temperature,  ammonia  and  trimethylgallium  will  rapidly  form  adduct  product: 
TMGa+NH3  4 ► TMGa:NH3  R-3 

The  remaining  of  the  mechanism  is  not  clear.  It  has  been  [54,55]  suggested  that  an  intra 
molecular  decomposition  path  is  followed: 


TMGa:NH3 


> DMGa:NH2  + CH4 


R-4 
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Then  this  monomer  can  form  a ring  structured  oligomer: 

xDMGa:NH2  ► (DMGa:NH2)x  R-5 

This  oligomer  may  then  by  an  unknown  mechanism  decomposes  on  the  substrate  surface 
to  form  GaN  film  [38]: 

(DMGa:NH2)x  ► GaN  product  R-6 

More  recently,  inter  molecular  reaction  mechanism  was  suggested  from  ab  initio 
calculations  [56],  In  this  mechanism,  adduct  TMGa:NH3  will  continue  react  with  TMGa 
then  NH3  to  produce  methane  to  form  a large  oligomer: 

TMGa:NH3  + TMGa  ► DMGa:NH2:TMGa  + CH4 

DMGa:NH2:TMGa  + NH3  ► (DMGa:NH2)2  + CH4  R-7 

Another  group  suggested  yet  another  possible  reaction  path  from  ab  initio 
calculation: 

TMGa:NH3  ► DMGa:NH2  + CH4 

DMGa:NH2  + TMGa  + NH3  ► (DMGa:NH2)2  + CH4 

R-8 

The  Sandia  National  Lab  GaN  group  reported  no  evidence  of  oligomer  formation 
based  on  their  Mass  spectroscopy  studies,  they  also  reached  the  conclusion  that  no 
TMGa:NH3  adduct  unimolecular  decomposition  occurs  [94].  But  they  did  not  rule  out 
high  order  reaction  mechanisms. 

The  gas  phase  chemistry  R-1  through  R-3  is  widely  accepted.  Experimental 
difficulty  have  limited  studies  on  the  subsequent  reactions  to  conducted  at  the  conditions 
far  away  from  the  real  growth  conditions.  Change  reaction  conditions  may  alter  the  gas 
phase  fluid  dynamics,  heat  and  mass  transfer  influences  on  the  gas  phase  chemistry. 
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In  this  contribution,  in  situ  Raman  spectroscopy  was  used  to  investigate  the  gas  phase 
chemistry  under  real  growth  conditions.  The  gas  phase  temperature  and  species 
concentration  as  a function  of  position  along  the  centerline  of  a cylindrical  inverted 
vertical  MOCVD  reactor  was  obtained  to  understand  the  gas  phase  reaction  mechanism. 


aZ 


Sweeping  flow  Reactant  flow 

Figure  4-1 . Inverted  impinging  flow  MOCVD  reactor. 

Experimental 

A 3-D  schematic  diagram  of  the  reactor  is  shown  in  Figure  4-1.  The  irmer  stream 
contains  the  Ga  source  (TMGa)  and  carrier  gas,  while  the  outer  stream  contains  the  NH3 
and  carrier  gas.  The  height  of  the  irmer  tube  can  be  adjusted  to  achieve  a different  degree 
of  mixing.  The  sweeping  stream  contains  carrier  gas  only.  In  the  TMGa+NH3+N2  study, 
V/III  =500  was  used.  For  the  TMGa;NH3+NH3+N2  study,  TMGa:NH3  was  obtained  by 
allowing  the  TMGa  and  NH3  to  flow  continuously  into  the  gas  distribution  zone,  which 
served  as  an  intermediate  container,  for  ten  hours  while  controlling  its  temperature  at  50 
°C.  The  gas  contribution  zone  was  packed  with  quartz  beads  of  3 mm  in  diameter  to 
provide  a in  line  container  for  condensation  of  the  adduct.  After  shutting  down  TMGa,  a 
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50%  mole  fraction  of  NH3  inN2  was  allowed  to  flow  continually.  The  adduct  was  then 
transported  under  its  own  vapor  pressure.  Data  were  taken  after  additional  30  min  gas 
flow  when  no  TMGa  signal  was  detected  at  the  gas  inlet.  Susceptor  temperature  was  kept 
at  1010  °C.  Raman  spectra  from  50  to  4500  cm‘‘  were  taken  along  the  centerline  of  the 
reactor.  Gas  phase  temperature  was  obtained  from  the  N2  rotational  spectra,  and  the 
species  concentration  profiles  were  calculated  as  described  in  Chapter  3,  from  the 
vibrational  transitions  when  Raman  cross-sections  were  available.  For  those  species  that 
Raman  cross-sections  are  not  available,  Raman  intensity  ratio  of  the  species  and  N2  were 


used  in  the  plot.  Standard  Raman  spectra  of  related  species  were  listed  in  Table  4-3. 


Table  4-3.  Vibrational  frequency  [cm'*]  of  known  species. 


Species 

C-Ga  Vs 

C-Nvs 

C-H  5s 

C-H  Vs 

N-H  Vs 

MMGa  [56] 

476 

1247.9 

2986.9 

— 

TMGa* 

526 

1210 

2919 

— 

TMGa:NHs  [57] 

544 

518 

1182 

2884 

3277 

(DMGa:NH2)3  [54]  540 

520 

1200 

2890 

3280 

CH4  [59] 

2917 

— 

C2H6[60] 

weak 

2954** 

— 

CH3  [61] 

3005 

— 

* Results  produced  in  this  lab. 

**  Usually  overlapped  with  C-H  Vas 

at  2968.69  cm 

-1 

Table  4-4.  ab  initio  calculation  results  [cm'*]  with  Gaussian98  [62] 

Species 

C-Ga  Vs  C-N  Vs 

C-H  5s 

N-H  5s 

C-H  Vs 

N-H  Vs 

MMGa 

472.3 

1246.4 

2986.9 

DMGa 

483.7 

1264.3 

3021.1 

TMGa 

511.6 

1297.1 

3019.6 

TMGa:NH3 

505.8  304.9 

1298.4 

1199.2 

2983.3 

3451.1 

MMGa:NH2 

495.8 

1257.2 

3003.6 

DMGa:NH2 

532.1  693.7 

1309.7 

1632.9 

3008.3 

3569.9 

(DMGa:NH2)2 

538.2  528.3 

1307.4 

1652.8 

2995.7 

3457.4 

ab  initio  calculation 


When  the  standard  spectrum  for  a species  has  not  been  published  or  produced  in 
our  own  laboratory,  ab  initio  calculation  results  must  used  as  an  alternative  source.  Table 
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4-4  lists  some  theoretical  spectra  obtained  from  ab  initio  calculations  using  Gaussian98 
[62],  For  comparison  the  ab  initio  calculation  results  of  MMGa,  TMGa  and  TMGarNHs 
are  also  listed. 

When  compare  the  ab  initio  calculation  results  with  experimental  results  that 
listed  in  the  Tables  4-3  and  4-4,  only  MMGa  has  an  excellent  match.  Farely  large 
discrepancy  was  found  for  TMGa.  Very  large  discrepancy  was  found  for  TMGa:NH3 
adduct.  Therefore,  the  ab  initio  calculation  results  at  the  current  stage  still  can  not  be 
used  for  large  molecule  species  identification. 

A 


10 


20 


I 

( 


> MMGa 



^ 

^ (DMGa:NH2)x 

A 


Flow 


30  1 


TMGa:NH3-NH3-N2 


TMGa:NH3 


Figure  4-2.  Identified  gas  phase  species  along  reactor  centerline. 

Results  and  Discussion 


TMGa:NH^+NHi+N9 

Figure  4-2  plots  the  identified  species  in  the  gas  stream  along  the  reactor 
centerline.  The  left  side  of  Figure  4-2  is  the  measured  gas  phase  temperature  profile  along 
the  center  line.  The  identification  of  these  species  was  straight  forward  except  that  of 
CH4.  The  C-H  stretching  line  in  the  range  2911  to  2919  cm’'  in  the  spectrum  of  TMGa  is 
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a superposition  of  the  vi(Ai’)  line  at  291 1 cm"'  and  the  vi2  (E)  at  2916  cm*'.  On  the  other 
hand  C-H  stretching  of  CH4  is  located  at  2917  cm"'.  Therefore,  it  is  difficult  to  separate 
CH4  from  TMGa.  The  identification  of  CH4  can  only  be  done  by  quantitatively 
calculating  of  the  intensity  of  TMGa  at  291 1-2919  cm’'  and  comparing  it  to  the  intensity 
of  TMGa  Ga-C  stretching  at  526  cm’’.  The  extraordinary  high  intensity  of  the  peak  at 
2911-2919  cm’'  suggests  of  the  existence  of  CH4. 
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Figure  4-3.  Gas  phase  species  concentration  profiles  along  the  reactor  centerline  for 
TMGa:NH3+NH3+N2  experiment. 


Carbon-Gallium  symmetric  stretching  (C-Ga  Vj)  after  deconvolution  were  used 
for  TMGa,  TMGa:NH3  and  (DMGa:NH2)x  concentration  calculations.  For  TMGa:NHx 
and  (DMGa:NH2)x  additional  information  from  Carbon-Nitrogen  symmetric  stretching 


Relative  Raman  Intensity 
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(C-N  Vs)  after  deconvolution  were  also  used.  Peak  positions  were  also  calibrated  against 
one  of  the  Ar^  plasma  line  at  530  cm''. 

Figure  4-3  shows  concentration  profiles,  where  the  TMGa  is  given  in  mole 
fraction  while  TMGa:NH3  and  oligomer  are  given  in  relative  Raman  intensity,  (ratio 
between  the  Raman  intensity  of  the  plotted  species  and  that  of  the  internal  reference  N2) 
since  the  Raman  cross-sections  of  these  species  are  unknown. 
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Figure  4-4.  Measured  gas  phase  temperature  profile  along  the  reactor  centerline. 

The  experimental  results  indicate  that  the  adduct  TMGaiNHa  begins  to  dissociate 
and  release  TMGa  when  the  gas  phase  temperature  begin  to  increase  at  the  edge  of  the 
thermal  boundary  layer.  The  TMGa  concentration  increases  rapidly,  peaking  at  around 
600  K,  then  decreases  as  the  probe  comes  closer  to  the  susceptor.  At  about  2 mm  beyond 
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the  location  of  first  detection  of  TMGa,  oligomer  becomes  detectable  and  increases 
monotonously  to  reaches  a plateau  when  temperature  is  greater  than  700K.  A possible 
explanation  of  this  behavior  is  that  the  oligomer  generation  reaction  rate  may  be  smaller 
than  TMGa  unimoleeular  decomposition  reaction  rate.  Figure  4-4  shows  the  measured 
gas  phase  temperature  profile  along  the  reactor  centerline  determined  from  the  N2 
rotational  Raman  transitions.  Figure  4-5  shows  the  measured  NH3  mole  fraction  profile 
along  the  reactor  centerline  determined  from  the  intensities  of  NH3  and  N2  vibrational 
Raman  transitions. 
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Figure  4-5.  Measured  NH3  mole  fraction  profile  along  the  reactor  centerline 
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Figure  4-6.  Measured  gas  phase  temperature  profile  along  the  reactor  centerline 
TMGa+NHci+N? 

Figure  4-6  shows  the  measured  gas  phase  temperature  profile  along  the  reactor 
centerline  determined  from  N2  rotational  Raman  transitions.  The  temperature  profile  is 
about  60  degree  lower  than  that  in  last  experiment  due  to  a layer  of  amorphous  material 
deposited  on  the  quartz  susceptor  that  reduced  the  thermal  radiation  to  the  inlet  and 
thermal  conduct  from  susceptor  to  the  gas  phase. 

Figure  4-7  shows  the  measured  NH3  mole  fraction  profile  along  the  centerline  as 
determined  from  the  intensities  of  NH3  and  N2  vibrational  Raman  transitions.  Figure  4-8 
shows  the  TMGa  concentration  profile  along  the  reactor  centerline  given  in  mole  fraction 
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while  the  TMGa:NH3  and  oligomer  profiles  are  given  in  relative  Raman  intensity.  Poor 
mixing  of  the  TMGa  and  NH3  allows  some  amount  of  TMGa  to  reach  the  detection  zone. 

Figure  4.8  shows  the  oligomer  in  the  low  temperature  region  increases  along  with 
TMGa  and  adduct  TMGa:NH3  decrease.  It  reaches  a plateau  when  TMGa  concentration 
drops  below  the  detection  limit.  The  concentration  increases  again  when  TMGa  reappear 
at  the  edge  of  the  thermal  boundary  layer  from  the  dissociation  of  adduct  TMGa;NH3.  Its 
concentration  reaches  yet  another  plateau  perhaps  due  to  the  same  mechanism  mentioned 
in  last  experiment.  The  lower  temperature  profile  results  in  a shorten  high  temperature 
region  (T>700K)  and  shorter  oligomer  high  concentration  plateau. 
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Figure  4-7.  Measured  NH3  mole  fraction  profile  along  the  reactor  centerline 
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Figure  4-8.  Gas  phase  species  concentration  profile  along  the  reactor  centerline 

Relevance  to  GaN  MOCVD 

From  these  two  experiments,  the  gas  phase  chemistry  likely  includes  the 
following  reactions: 


NH3  -I-  M 


NH2  + H 


G-1 


ki  k2  k3 

TMGa  ► DMGa  ► MMGa  ► Ga  G-2 
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TMGa:NH3 
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G-3 

DMGa:NH2:TMGa  + CH4 
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DMGa:NH2:TMGa  + NH3  ► (DMGa:NH2)2  + CH4 

G-4 

From  the  oligomer  concentration  profiles  showed  in  Figure  4-3  and  Figure  4-8, 
the  two  regions  where  oligomer  concentration  increases  rapidly  are  appearent:  the  inlet 
region  and  the  edge  of  the  thermal  boundary  layer  region.  Reactor  design  and  operation 
conditions  should  avoid  large  amounts  of  oligomer  formation  in  these  two  regions. 

At  inlet  region,  reduce  reaction  G-4  by  allowing  no  adduct  formation  reaction  or 
allowing  adduct  formation  reaction  to  complete  by  enhancing  the  forward  reaction  in 
G-3.  This  can  be  done  by  using 

• very  high  V/III  ratio, 

• high  pressure, 

• low  temperature, 

• long  residence  time. 

Since  adduct  TMGaiNHs  has  very  low  vapor  pressure,  reactor  must  be  operated  in 
super  saturation  condition.  The  reactor  wall  must  be  warmed  up  to  allow  condensed 
adduct  TMGa:NH3  evaporate  back  to  the  gas  phase  while  the  wall  temperature  must  be 
lower  enough  to  avoid  its  dissociation. 

In  the  high  temperature  region,  the  delivery  of  small  molecules  to  the  surface  is 
promoted  by  creating  steep  increasing  temperature  profile,  promotes  the  reverse  reaction 
in  G-3,  reactions  G-1  and  G-2  to  occur  inside  the  thermal  boundary  layer  where  the  gas 
phase  temperature  is  higher  than  700  K. 

Device  quality  GaN  is  usually  grown  at  high  substrate  temperature  of  1050  °C. 
With  a steep  temperature  gradient,  buoyant  force  can  be  important  role  in  an  impinging 
flow  reactor.  The  Richardson  number,  the  relative  importance  of  buoyant  forces  to 
forced  convection  is  defined  as  follows. 
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Ri  = (Ap  / p)gL  / ul  4-1 

At  high  Ri  condition  (>1),  the  buoyancy  effect  will  dramatically  alter  the 
temperature  profile.  Therefore,  to  reduce  the  buoyancy  effect,  high  flow  velocity  or  low 
concentration  of  reactants  are  needed,  which  will  be  discussed  in  detail  in  Chapter  7. 
From  Equation  4- 1 , the  steep  temperature  profile  in  impinging  flow  can  be  realized  by, 

• high  velocity, 

• N2  as  carrier  gas, 

• two  flow  design, 

• high  speed  rotating  disk  design. 

These  conditions  are  very  familiar  since  many  of  them  have  already  been  incorporated 
into  the  reactor  design  and  MOCVD  operation  through  trial  and  error. 

The  conclutions  of  this  study  are: 

• The  gas  phase  temperature  profile  important  in  determine  the  reaction  path. 

• An  alternative  gas  phase  reaction  path  suggested. 

• The  adduct  oligomer  (DMGa:NH2)x  formed  by  reaction  of  the  parent  adduct 
TMGa:NH3  with  TMGa  and  NH3  is  revealed  in  the  Raman  experiments. 


CHAPTER  5 

PARAMETER  IDENTIFICATION  IN  MOCVD: 

USING  IN  SITU  RAMAN  SPECTROSCOPY, 

CVD  MATHEMATICAL  SIMULATION 
AND  INVERSE  PROBLEM  FORMULATION 

Introduction 

Mathematical  CVD  models  built  upon  the  fundamental  transport  phenomena 
governing  equations  have  been  used  by  many  research  groups  [7-1 1,63-66]  to  provide 
accurate  simulation  of  fluid  patterns,  temperature  and  species  concentration  profiles  with 
complex  chemical  reactions.  The  success  was  dictated  by  the  ability  to  obtain  a clear  and 
reliable  understanding  of  the  reaction  mechanisms  as  well  as  chemical  and  physical 
properties  of  all  related  species.  Experimental  data  for  many  of  these  properties  may  not 
be  available  especially  for  those  newly  emerged  compound  semiconductor  materials. 
Often  time  theoretical  estimation  was  used.  On  the  other  hand  spectroscopic 
measurements  can  provide  a wide  range  of  on  temperature  and  species  concentrations  in 
real  reaction  conditions.  Thus  it  is  natural  to  link  measurements  and  modeling  together  to 
perform  an  inverse  simulation  from  which  species  chemical  and  physical  properties  may 
be  identified.  With  the  property  values,  process  optimization  and  reactor  design  may 
then  be  performed. 

In  this  chapter  a methodology  will  be  extablished  to  combine  in  situ  Raman 
spectroscopy  with  numerical  CVD  modeling  and  inverse  problem,  from  which  species 
chemical  and  physical  properties  can  be  identified.  Uniqueness  of  the  problem  will  also 
be  examined,  and  a numerical  solution  method  will  be  tested. 
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Inverse  Problem:  Theoretical  Framework 
A discrete  parameter  identification  problem  takes  the  general  form  of 

min||w(a)-2‘'||^ . 5_1 

For  a more  specific  case  such  as  elliptic  second  order  differential  equation,  constraints 
may  apply 

Au-  f in  Q 5-2 

with  the  Dirichlet  boundary  condition 

u = g on  5Q  5-3 

and  additional  boundary  data 

du  / dN ^ =h  on  r e 5Q  5-4 

The  problem  is  often  referred  to  as  an  ill-posed  problem  [18].  Therefore,  before  solving 
the  problem,  the  illness  issues  must  be  tackled  first. 

Uniqueness  of  the  problem  using  single  boundary  measurement  data 

The  uniqueness  of  identifying  the  parameters  of  an  elliptic  partial  differential 
equation  was  mathematically  proved  by  Robert  Kohn  and  Michael  Vogelius  [67].  It 
states  that  these  parameters  can  be  uniquely  determined  by  the  so-called  Dirichlet-to- 
Neumarm  mapping  that  allows  only  the  measurements  along  the  centerline  of  a 
cylindrical  domain  (called  bore  measurements).  Essentially,  the  Dirichlet-to-Neumann 
mapping  calls  for  each  iteration  step,  one  simulation  with  measured  concentration  profile 
then  another  simulation  with  measured  flux  profile  on  the  centerline,  then  uses  the 
calculated  differences  of  the  simulations  to  update  the  identifying  parameter.  This  is  a 
very  robust  method  but  obviously  costly. 
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Another  important  contribution  was  that  by  A.  Khaidarov  [68],  He  has  proved 
that  the  parameters  of  an  elliptic  equation  can  be  uniquely  determined  by  any  single 
boundary  measurements  if  these  parameters  are  independent  of  at  least  one  coordinate 
variable.  His  proof  indicates  that  measurements  on  a single  boundary  will  contain 
enough  information  for  a unique  recovery  of  all  parameters,  this  makes  it  possible  to  use 
a much  less  costly  method,  the  so-called  Output  Least  Square  (OLS)  method  with 
Tikhonov  regularization. 

Output  Least  Square  (OLS)  Method  with  Tikhonov  Regularization 

The  noisy  nature  of  the  experimental  observations  will  destabilize  the  system  thus 
causing  oscillation.  The  Tikhonov  regularization  method  was  originally  suggested  by 
Tikhonov  for  dealing  with  the  unstable  nature  of  linear  ill-posed  problems  [18],  and 
recently  it  have  been  adopted  for  solving  nonlinear  problems  [69,70].  Since  the  Navier- 
Stokes  equation,  energy  conservation  equation  and  species  conservation  equations  are 
usually  categorized  as  weak  nonlinear,  Tikhonov  regularization  method  should  be 
adequate  for  solving  this  problem. 

In  Tikhonov  regularization  method.  Equation  5-1  is  replaced  by 


where  x and  y represent  the  Hilbert  spaces;  u{a)  is  the  mole  fraction  from  a simulation 
and  a function  of  reaction  rate  constant  a\  z®  is  the  measured  mole  fraction;  (5  is  the 
Tikhonov  regularization  parameter;  P represent  the  projection  from  Hilbert  space  x to 
Hilbert  space  y’.  This  projection  is  designed  to  allow  the  parameters  approach  to  its  true 
value  at  the  same  speed  as  its  counterpart  dose  to  reduce  the  residue  error. 


5-5 
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Figure  5-1  A schematic  flow  chart  for  the  parameter  identification  algorithm. 

The  schematic  flow  chart  for  the  parameter  identification  algorithm  is  shown  in 
Figure  5-1.  Equation  5-5  is  equivalent  to  the  Levenberg-Marquardt  method  [71] 

(j’A  + A^/)x^  ^A’b  5-6 


where  / is  the  identity  matrix. 

Reactor  Simulation  — the  Forward  Problem 

The  reactor  model  is  axisymmetric,  steady  state,  and  consists  of  momentum, 
energy  and  mass  conservation  equations  in  cylindrical  coordinates  [5]. 

R-direction  momemtum  conservation  equation: 


dv^ 
P(vr  ^ 
or 


dvy.  ^ dp 

+ Vz  — ^)  = — - + P 
^dz  dr 


dr 


1 d 

(rv  ) 

rdr  ^ , 


d V 


+ • 
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Z-direction  momemtum  conservation  equation: 


^ dvz  dvz^  dp 

Piyr  ^ + Vz-— ) = - — + // 
dr  dz  dz 


1 d ^ 
r dr  dr 


d V 
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Energy  conservation  equation: 


dT 


dT. 


pc  ivy.  — + vz  —)^k 


dr 


dz 


1 d dT  d^T 
{r — ) + — ^ 


r dr  dr 


dz‘ 
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Species  conservation  equation: 


Sx. 


dx. 


1 d 


P(yr  ^ + vz  -^)  = , 

dr  dz  r dr 
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+ =0 


Mass  conservation  equation: 


r dr  ' dz  ^ 


5-11 


The  boundary  conditions  of  momentum  and  energy  conservation  equations  are  shown  in 
Figures  5-2  and  5-3.  Finite  element  method  was  used  in  the  forward  CVD  modeling, 
while  the  Levenberg-Marquardt  [70]  method  was  used  for  updating  the  parameter. 

NHt  Homogeneous  Decomposition  Rate  Constant  Identification 

Introduction 

NHa-Ar-system  has  been  studied  by  several  groups  using  so-called  shock  tube 
method,  in  which  a second  order  reaction  has  been  identified  and  the  first  step  of  the  gas 
phase  reaction  was  suggested  as  follow: 


k, 

NH3  + Ar  > NH2  + H + Ar  R5-1 

The  reported  reaction  rate  constants  00  the  first  step  ki  by  different  groups  all 


summarized  in  Table  4-2. 
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All  reactor  inner  walls: 
Non-slip  condition 


Figure  5-2.  Schematic  momentum  equation  boundary  conditions  of  the  reactor  model 


Figure  5-3.  Schematic  energy  equation  boundary  conditions  of  the  reactor  model. 
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NHa  heterogeneous  reaction  on  quartz  surface  was  analyzed  by  Dogh  Myer  [72], 
Reaction  mechanism  was  suggested  as  follow: 

NH3  + S > NH-S  + H2  +S  R5-2 

Overall  reaction  rate  constant  was  also  identified. 


Heater: 

Constant  temperature 


On  substrate: 

Flux  equal  to  the  rate  of 
heterogeneous  Reaction 

Inlet: 

Constant  concentration 


All  reactor  inner  walls: 
Species  flux  equals  zero 


Figure  5-4.  Schematic  concentration  boundary  conditions  of  the  NH3  decomposition 

Spatially  resolved  gas  phase  temperature  and  species  concentration  profiles  along 
the  centerline  of  a vertical  MOCVD  reactor  were  obtained  by  rotational  Raman 
temperature  measurement  and  vibrational  Raman  species  concentration  measurement 
technique  described  in  Chapter  3.  The  reasonable  minimum  signal  to  noise  ratio  S/N 
commonly  accepted  for  a quantitative  interpretation  is  about  2.0.  The  S/N  for  any 
intermediate  or  product  in  both  cases  was  not  large  enough  to  have  a meaningful 
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quantitative  interpretation.  Therefore,  the  gas  phase  reaction  rate  constant  identification 
will  only  limited  to  the  gas  phase  reactant  disappearing  rate  at  this  stage. 

Experimental 

The  gas  phase  reaction  system  with  in  situ  Raman  spectroscopy  described  in 
Chapter  11  was  used  for  this  study.  High  purity  N2  (99.999%)  was  introduced  into  the 
center  tube.  50%  volume  ratio  of  blue  NH3  (99.999%)  in  N2  was  introduced  into  the 
annular  region  of  the  center  tube  packed  with  3mm  diameter  quartz  beads.  Pure  N2  was 
used  as  sweeping  flow.  All  inlets  flow  velocity  was  set  at  2.0  [cm/s].  Reactor  aspect 
ratio  was  set  at  0.5,  and  the  flow  velocity  was  calculated  at  296.15  K and  one 
atmospheric  pressure. 


Figure  5-5  Center  line  gas  phase  temperature  profile. 
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The  heater  temperature  was  set  at  1373.15  K measured  by  the  online 
thermocouple.  The  entire  reactor  was  kept  at  one  atmospheric  pressure.  Raman  spectra 
were  taken  along  the  center  line  of  the  reactor.  The  N2  rotational  part  of  the  Raman 
spectra  was  used  for  gas  phase  temperature  measurement,  while  for  N2  at  233 1 cm''  and 
for  NH3  at  3337  cm''  vibrational  Raman  transition  peaks  were  used  for  the  NH3 
concentration  profile  measurement.  Details  of  the  data  analysis  procedures  are  described 
in  Chapter  3.  The  NH3  heterogeneous  reaction  at  quartz  surface  proposed  by  Meyer  [72] 
was  used. 


Figure  5-6  Center  line  ammonia  mole  fraction  profile 
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Results  and  Discussion 

Figure  5-6  shows  the  experimental  and  simulated  gas  phase  NH3  concentration 
profiles  along  the  center  line.  The  simulation  also  includes  two  hypothetical  cases:  1)  no 
reaction,  2)  gas  phase  reaction  only  that  often  suggested  in  the  literatures.  The  gas  phase 
temperature  profile  is  shown  in  Figure  5-5. 


Temperature  [K] 
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Figure  5-7.  Simulated  flow  steam  lines  and  the  spatial  distribution  of  gas  phase 
temperature  and  NH3  concentration  of  the  entire  flow  region 


The  umerical  procedures  described  in  the  earlier  sections  of  this  chapter  were 
used  for  NH3  gas  phase  thermal  decomposition  rate  constant  identification  using  the 
experimental  results  obtained  by  in  situ  Raman  spectroscopy,  under  the  assumption  of 
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second  order  homogeneous  reaction.  As  shown  in  Figures  5-5  and  5-6,  the  simulation 
results  of  both  gas  phase  temperature  and  the  NH3  concentration  profile  along  the  center 
line  of  the  reactor  fit  the  in  situ  Raman  experimental  results  very  well.  The  reaction  rate 
constant  of  the  first  step  of  NH3  thermal  decomposition,  ki,  is  therefore  was  found  to  be 
1.4xl0'^exp(90.1 1 kcal/RT)  [cm^  mof'  s’'].  The  identified  value  matches  well  with  the 
literature  reported  values  listed  in  Table  4-2. 

From  left  to  right  side  the  simulated  flow  steam  line,  and  spatial  distribution  of 
gas  phase  temperature  and  NH3  concentration  of  the  entire  flow  region  are  shown  in 
Figure  5-7. 


TMGa-N?  Binary  Diffusion  Coefficient  Identification 

Introduction 

Experimental  measurement  of  TMGa-N2  binary  diffusion  coefficient  has  never 
been  reported.  The  binary  diffusion  coefficient  can  be  estimated  using  kinetic  theory  [5]. 
The  error  by  this  method  is  often  within  3 to  5%  [5].  The  mean  collision  diameter  cr  and 
sk  of  TMGa  were  reported  [66]  to  be  5.470  A and  378.2  K,  respectively.  Heat  capacity 
Cv  was  estimated  in  this  work  to  be  26.271  [cal/mol-K]  calculated  using  Gaussian  98 
[62]. 

Binary  diffusion  coefficient  in  the  order  Chapman-Enskog  approximation  takes 
the  form  [74], 


[^2  ]i 


3 

\6nm^2 


5-12 


Therefore,  the  theoretical  estimation  of  this  binary  mass  diffusion  coefficient,  [S>i2]i,  is 
2 1 

0.08024  [cm  s’  ].  The  usually  3 to  5%  errors  can  be  improved  through  application  of 
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Chapman-Cowling  approximation.  The  2"**  order  Chapman-Cowling  approximation  for 
the  binary  diffusion  coefficient  can  be  written  in  the  form  [74] 

[%]2=[^2],  -/d  5-13 

where /o  =(1-A,2)“'  5-14 
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5-23 

8^2 

From  the  equations  above,  one  can  find  that  the  3 to  5%  errors  are  concentration 
dependent.  This  concentration  dependency  will  also  undermine  the  uniqueness  of  the 
identification  process  therefore  must  be  removed.  By  assuming  the  true  binary  diffusion 
coefficient  equals  to  the  order  Chapman-Cowling  approximation  of  the  binary 
diffusion  coefficient,  Equation  5-13  can  be  used  to  deconvolute  the  binary  diffusion 
coefficient  allow  the  identification  of  the  concentration  independent  order  Chapman- 

Enskog  approximation. 

Experimental 

The  gas  phase  reaction  system  with  in  situ  Raman  spectroscopy  described  in 
Chapter  2 was  used  for  this  study.  High  purity  N2  (99.999%)  carrier  gas  at  flow  velocity 
of  2.0  cm/s  was  introduced  into  the  center  tube,  along  with  5%  mole  fraction  of  high 
purity  TMGa  (99.999%)  as  a reactant.  Pure  N2  at  the  same  flow  velocity  was  introduced 
into  annular  region  of  the  center  tube  and  the  sweeping  flow  regions.  Both  were  packed 
with  3mm  diameter  glass  beads.  Reactor  aspect  ratio  was  set  at  0.5.  Flow  velocity  was 
calculated  at  296.15  K and  one  atmospheric  pressure.  The  heater  temperature  was  set  at 
room  temperature  and  the  entire  reactor  was  kept  at  one  atmospheric  pressure.  Raman 
spectra  were  taken  along  the  center  line  of  the  reactor.  N2  at  233 1 [cm’']  and  TMGa  at 
526  [cm"']  vibrational  Raman  transition  peaks  were  used  for  TMGa  concentration  profile 
measurement.  Details  of  the  data  process  procedures  are  described  in  Chapter  3.  By 
following  the  description  the  TMGa  binary  diffusivity  are  also  being  estimated. 
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Results  and  Discussion 

The  numerical  procedures  described  in  the  earlier  sections  of  this  Chapter  were 
used  to  identify  TMGa-N2  binary  mass  diffusion  coefficient  using  the  experimental 
results  obtained  by  in  situ  Raman  spectroscopy.  As  shown  in  the  Figure  5-8,  the 
simulation  results  of  TMGa  concentration  profile  at  the  center  line  of  the  reactor  fit  the  in 
situ  Raman  experimental  results  very  well.  The  binary  mass  diffusion  coefficient  was 
found  to  be  as  [S>i2]i-0.08  [ernes''],  which  agrees  very  well  with  the  theoretical  estimated 
value  0.08024  [ernes'*].  The  order  Chapman-Cowling  approximation  correction  term 
is  decovoluted,  not  include  in  this  value. 


Figure  5-8.  Reactor  Center  Line  TMGa  concentration  profile. 
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TMGa  [w%] 


Figure  5-9.  Simulated  flow  steam  line  and  spatial  distribution  of  gas  phase  TMGa 
concentration 

As  shown  in  Figure  5-8,  TMGa  mole  fraction  drops  from  0.05  to  about  0.002  within  6 
[mm]  from  the  inlet.  This  rapid  concentration  change  is  not  due  to  the  diffusion  since  the 
diffusion  coefficient  is  very  small  compare  with  ammonia.  The  true  reason  for  this  rapid 
concentration  drop  is  due  to  buoyancy  effect.  Because  the  high  molecular  weight  of 
TMGa,  the  density  gradient  is  sufficient  for  make  the  flow  unstable.  Predictably,  if 
hydrogen  was  used  as  the  carrier  gas,  the  buoyancy  effect  would  be  even  more  dramatic 
since  the  TMGa  diffusion  coefficient  and  molecular  weight  difference  are  greater.  Figure 
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5-9  gives  the  better  visual  of  this  phenomenon;  TMGa  convected  transverse  to  the  flow 
axis. 


Heater: 

Constant  temperature 


4 


Inlet: 

Constant  concentration 


Figure  5-10.  Schematic  concentration  boundary  conditions  of  the  TMGa  decomposition 
reaction  rate  constant  identification 

TMGa-N?  Thermal  Decomposition 

Introduction 

This  reaction  has  been  studied  extensively,  by  both  ab  initio  quantum  mechanical 
calculations  and  experimental  measurements  have  been  reported  by  several  researchers. 
A gas  phase  decomposition  mechanism  was  also  proposed  [42-45]  and  reaction  rate 
constants  for  each  step  have  been  reported  [42,43].  It  is  clear  that  trimethylgallium 
(TMGa)  decompose  to  form  dimethylgallium  (DMGa),  then  monomethylgallium 
(MMGa): 
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ki  lc2  ks 

TMGa  DMGa  ^ MMGa  ^ Ga  R-1 

The  reported  rate  constants  are  summarized  in  Table  4-1 . 


Figure  5-11.  Reactor  center  line  phase  temperature  profile 
Experimental 

Gas  phase  reaction  system  with  in  situ  Raman  spectroscopy  described  in 
Chapter  2 was  used  for  this  study.  High  purity  N2  (99.999%)  at  a flow  velocity  of  3.0 
cm/s  was  introduced  into  the  center  tube  packed  with  3mm  diameter  quartz  beads  as 
shown  in  Figure  5-10.  This  stream  contains  5%  weight  fraction  of  high  purity  TMGa 
(99.999%)  in  N2  carrier  gas  as  reactants  high  purity  N2  (99.999%)  was  used  as  sweeping 
flow.  Reactor  aspect  ratio  was  set  at  1 .0.  Flow  velocity  was  calculated  at  296. 1 5 K and 
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one  atmospheric  pressure.  The  heater  temperature  was  kept  at  1373.15  K measured  by  an 
online  thermocouple  and  the  entire  reactor  was  kept  at  one  atmospheric  pressure.  Raman 
spectra  were  taken  along  the  center  line  of  the  reactor.  N2  rotational  part  of  the  Raman 
spectra  was  used  for  gas  phase  temperature  measurement  while  N2  at  233 1 [cm"']  and 
TMGa  at  526  [cm'']  vibrational  Raman  transition  peaks  were  used  in  these 
measurements.  Details  of  the  data  reduction  procedure  were  described  in  Chapter  3. 


Figure  5-12.  Reactor  center  line  TMGa  concentration  profiles 
Results  and  Discussion 

Figure  5.1 1 shows  the  experimental  and  simulated  gas  phase  temperature  profiles 
along  the  center  line  the  reactor.  The  TMGa  concentration,  as  shown  in  the  Figure  5-12, 
decreases  rapidly  in  the  high  temperature  region.  Simulation  results  of  both  gas  phase 
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temperature  and  TMGa  concentration  profiles  fit  in  situ  Raman  experimental  results  very 
well.  The  TMGa  decomposition  reaction  rate  constant  ki  was  found  to  be  3.4xl0'^ 
exp(60.1  kCal/RT)  [s  ‘],  which  matches  very  well  with  the  literature  value  [42]  listed  in 
Table  4-1.  The  simulated  flow  steam  line,  and  spatial  distribution  of  gas  phase 
temperature  and  TMGa  concentration  of  the  entire  flow  region  are  shown  in  Figure  5-13. 


Temperature  [K] 


TMGa  [w%] 


1223.46 

1014.56 

805.65 

596.744 

387.837 


0.146129 

0.116272 

0.0864147 

0.0565576 

0.0267006 


Figure  5-13.  Simulated  flow  steam  line,  and  spatial  distribution  of  gas  phase  temperature 
and  TMGa  concentration  of  the  entire  flow  region 


CHAPTER  6 

MULTICOMPONENT  DIFFUSION 
AND  THERMAL  DIFFUSION  COEFFICIENT 

Almost  all  CVD  simulations  involve  the  modeling  of  equation  of  changes  for  a 

multicomponent  system.  In  a multicomponent  system,  according  to  irreversible 

processes,  there  will  be  a contribution  to  each  flux  owing  to  each  driving  force  in  the 

system,  therefore  coupling  may  occur  [5].  Consequently,  in  a multicomponent  system, 

the  momentum  flux  depends  only  upon  the  velocity  gradients,  the  energy  flux  depends 

both  on  the  temperature  gradient  (heat  conduction)  and  on  the  mechanical  driving  forces 

(the  “diffusion-thermo  effect”  or  “Dufour  effect”),  and  the  mass  flux  depends  both  on  the 

mechanical  driving  forces  (concentration  gradient,  pressure  and  forced  diffusion)  and  on 

the  temperature  gradient  (the  “thermal-diffusion  effect”  or  “Soret  effect”).  The  Onsager 

reciprocal  relations  of  the  thermodynamics  of  irreversible  processes  give  information  as 

to  the  interrelation  of  the  two  coupled  effects.  Therefore  in  addition  to  transport 

properties  such  as  viscosity,  thermal  conductivity  and  diffusivity,  the  “thermal-diffusion 

ratio,  kj”  needs  to  be  introduced  to  describe  the  thermal  diffusion  effect.  Because  of  the 

interconnection  of  the  Soret  and  Dufour  effects,  as  described  by  the  Onsager  relations, 

the  “thermal-diffusion  ratio,  ^t,”  will  be  able  to  take  care  of  both  phenomena  [5], 

Thermal  Diffusion  and  Diffusion  Thermo 

Thermal  diffusion  is  caused  by  the  relative  motion  of  the  components  of  a mixture 

due  to  the  presence  of  a temperature  gradient.  As  a result  of  this  motion,  composition 

gradients  subsequently  appear  in  the  mixture  that  produces  ordinary  diffusion  that  tends 
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to  eliminate  these  gradients.  A steady  state  may  be  finally  reached  in  which  the 
separation  effect  arising  from  thermal  diffusion  is  balanced  by  the  remixing  effect  of 
ordinary  diffusion.  Partial  separation  is  then  observed  with  the  heavy  components 
usually  in  the  colder  region  and  the  light  components  in  the  hotter  region.  This 
phenomenon  was  first  observed  in  liquids  by  Ludwig  and  later  by  Soret,  and  is  now 
known  as  the  Soret  effect.  In  the  case  of  gases,  although  the  existence  of  thermal 
diffusion  had  been  suspected  by  Fedderson  already  in  1873,  it  was  predicted  theoretically 
by  Chapman  and  Enskog  almost  simultaneously  and  later  confirmed  in  the  experiments 
by  Chapman  and  Dootson. 

In  general,  the  thermal-diffusion  and  the  diffusion-thermo  effects  are  of  a smaller 
order  of  magnitude  than  the  effects  described  by  Fourier’s  or  Fick’s  laws  and  are  often 
neglected  in  the  calculation  of  heat  and  mass-transfer  process.  There  are,  however, 
exceptions.  In  a typical  GaN  MOCVD,  as  discussed  in  Chapter  4,  in  order  to  control  the 
desired  gas  phase  reaction  path,  a very  large  temperature  gradient  and  very  low  metal 
organic  species  concentration  are  used.  There  also  exist  a very  large  molecular  weight 
difference  between  carrier  gas  (H2  or  N2)  and  metal  organic  species  (e.g.  TMGa,  adduct 
and  oligomer).  These  special  conditions  can  largely  enhance  the  Soret  effect  to  cause  it 
to  be  of  a magnitude  such  that  it  cannot  be  neglected  [5]. 

History  of  Diffusion  Studies 

Fick’s  Law 

1855  Fick  proposed  a diffusion  law,  usually  known  as  “Fick’s  first  law”  for 
binary  mixture. 


i 
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where  ^ab  is  called  the  Pick’s  law  binary  diffusion  coefficient.  However  there  was  no 

precise  statement  as  to  which  flux  (left  hand)  was  to  be  set  to  equal  to  (right 

hand)  until  Bird  [5]. 

Maxwell  and  Stefan 

Developed  independently  by  Maxwell  and  Stefan,  a theory  of  momentum 
transport  that  led  to  an  expression  for  diffusive  fluxes  in  a isobaric  binary  system  in  the 
form  of  a partial  momentum  balance  (later  found  to  be  also  good  for  nonisobaric 
isothermal  binary  mixture), 

_ Fa^b  -Fb°^a 
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where  ^ab  is  called  the  Maxwell-Stefan  diffusivity. 

They  also  predicted  that  the  binary  molecular  diffusivity. 
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Generalized  Maxwell-Stefan  Formulation 

For  non-ideal  mixtures  in  the  presence  of  simultaneous  heat  transfer,  a 
generalized  Maxwell-Stefan  (GMS)  formulation  was  introduced.  The  relationship  of 
GMS  diffusivity  and  Pick’s  Law  diffusivity  was  found  as. 


2)2  =^2 
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where  yi  is  the  activity  coefficient  in  solution,  therefore,  for  ideal  gas  mixtures,  the  GMS 
diffusivity  will  reduce  to  Pick’s  law  diffusivity  • This  allows  the  application  of 

the  interpretation  of  diffusivity  by  molecular  collision  theory. 
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Thermodynamics  of  Irreversible  Processes 

The  fundamental  postulate  of  the  thermodynamics  of  irreversible  processes  states 
that,  for  systems  which  are  not  too  far  from  equilibrium,  the  fluxes  can  be  expressed  as 
r.  6-5 

j 

where  the  coefficients  «jj  are  phenomenological.  This  equation  states  that  the  fluxes  are 
linear  functions  of  the  driving  forces.  The  fundamental  theorem  of  the  thermodynamics 
of  irreversible  processes,  due  to  Onsager  states  that,  for  a suitable  selection  of  the  T j and 
X\,  the  phenomenological  coefficients  are  symmetric: 

a,j  = . 6-6 

This  is  known  as  Onsager  reciprocal  relations. 

Multicomponent  Species  Equations  of  Change 
Among  the  many  published  mathematical  CVD  simulation  modeling,  some  of 
them  use  commercial  computational  fluid  dynamics  (CFD)  codes  that  use  the  so  called 
effective  binary  approximation  method,  which  was  found  to  have  up  to  70%  errors  [76]. 
Some  of  them  use  more  elaborate  formulations  that  suppose  to  allow  more  accurate 
results.  However,  recent  progress  made  in  nonequilibrium  thermodynamics  theory  for 
coupled  heat  and  mass  transport  [77,78]  led  to  Curtiss  and  Bird  to  reexam  this  subject 
[79].  They  pointed  out  that  the  molecular  theory  led  to  expressions  for  the  mass  fluxes  of 
the  Maxwell-Stefan  form  but  also  predicted  that  the  mass  fluxes  would  be  dependent  on 
the  velocity  gradients  in  the  system.  Such  dependence  is  not  allowed  in  classical 
irreversible  thermodynamics,  where  only  linear  flux-force  relations  are  studied.  This  has 
led  them  to  conclude  that  some  earlier  as  well  as  current  work  (following  the  references 
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therein)  were  not  consistent  with  the  thermodynamics  of  irreversible  processes  for  linear 
systems  [69]. 

Since  their  early  formulation  [5,73]  formed  the  backbone  of  most  transport 
phenomena  simulation  especially  CVD  modeling,  CVD  models  had  go  also  be 
reformulated.  In  this  chapter,  following  this  reformulation,  a new  multicomponent  CVD 
model  was  developed.  This  model  may  be  used  for  future  parameter  identification  and 
process  optimization  in  multicomponent  system. 

The  species  equations  of  change  of  a multicomponent  system  is  as  [5] 
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and 
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where 


V-p,y  = Vpi  -Y  + yO^CV-y). 
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From  the  incompressible  fluid  assumption  V- v = 0 [5],  therefore. 


— n - -Vp.  • V - V • j + r 
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From  to  Curtiss  and  Bird’s  formulation  [79], 
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This  can  be  simplified  to  have  the  form: 
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and 
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From  Ferziger  [74]  the  multicomponent  thermal  diffusion  coefficients  are  given  as, 
Df  =p,Df  =p.^Dykj  6-14 

;=i 

therefore, 
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Multicomponent  Diffusion  and  Thermal  Diffusion  Coefficients 
Multicomponent  Diffusion  Coefficients 

From  Curtiss  and  Bird’s  new  formulation  [79],  the  multicomponent  diffusion 
coefficients  can  be  calculated  from 

ZZ[{^«*(l-<^,*)-<^,*  ](Oj -5„jCO,]c,„D,^j  ^Q),cOj  j ^ i 6-18 

k^j  nH 


where 
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which  has  the  symmetric  property 

Cy=Cj,,  6-20 

where  is  the  Pick’s  law  binary  diffusion  coefficient  [79],  Therefore  the 
multicomponent  diffusion  coefficients  have  the  properties 
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which  has  the  symmetric  property, 
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Similar  to  that  in  the  Chapter  5,  the  binary  diffusion  coefficient  in  the  1^'  order 
Chapman-Enskog  approximation  takes  the  form  [74], 
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In  the  2"^*  order  Chapman-Cowling  approximation  the  binary  diffusion  coefficient  can  be 
written  in  the  form  [74] 


i^j]2  ~ [^y]|  ■ fo 
where  for  multicomponent 

/o=(i-[A,Lr' 

and 
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and 
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Again  the  coefficient  A”  must  use  the  formulation  provided  by  Ferziger  [74]  to  guarantee 
the  symmetric  nature  of  the  diffusion  coefficients 

• 6-28 
Multicomponent  Thermal  Diffusion  Coefficients 


From  Ferziger  [74],  the  multicomponent  thermal  diffusion  ratios  are  as  follow 


and 
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where  x - M — - is  the  molar  fraction  of  species  / and  cOj  is  the  mass  fraction  of  species 


i.  Coefficient  A”  is  the  partial  bracket  integrals  can  be  written  as  a linear  combination  of 
the  Q-integrals,  the  diffusion  functional,  given  in  Ferziger  [74],  Here  it  must  be  point  out 
that  coefficient  A”  must  use  the  formulation  provided  by  Ferziger  not  the  traditional 

formulation  provided  by  Hirschfelder,  Curtiss  and  Bird  [73]  to  have  the  desired 
symmetric  character.  Therefore  the  thermal  diffusion  coefficients  will  have  the 
properties 
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J 

The  equations  deduced  above  completely  define  the  multicomponent  diffusion  and 
thermal  diffusion  coefficients  for  the  multicomponent  species  conservation  equation  used 
in  this  study,  which  incorporated  the  most  current  achievements  in  this  area.  And  the 
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multicomponent  diffusion  coefficients  are  symmetric.  Species  weight  fraction  used  for 
calculation  in  Table  6-2  is  given  in  Table  6-3. 


Table  6-1  Lennard-Jones  parameters 


species 

MW(g/mole) 

cr(A) 

(K) 

Ref 

H2 

2.0158 

2.968 

33.3 

[80] 

NHs 

17.0317 

2.915 

220 

[81] 

TMGa 

114.8211 

3.4 

378.2 

[82] 

TMGaiNHs 

131.8528 

5.47 

278.9 

[82] 

(DMGa:NH2)3 

347.4336 

5.89 

273.2 

[82] 

N2 

28.016 

28.016 

98.4 

[80] 

2 

Table  6-2.  Calculated  diffusion  and  thermal  diffusion  coefficients  [cm  /s] 


i-j 

[^]i 

[^]2 

[4]2 

D, 

XlO^ 

kij 

H2-H2 

0.0000 

0.0000 

0.0243 

-2.72x10'^ 

H2-NH3 

0.7269 

0.7394 

0.0169 

-5.827 

H2-TMGa 

0.3419 

0.3432 

0.0039 

-5.082 

H2-TMGa:NH3 

0.3449 

0.3499 

0.0143 

-5.098 

H2-(DMGa:NH2)3 

0.2408 

0.2445 

0.0152 

-4.921 

H2-N2 

0.7550 

0.7777 

0.0292 

-5.894 

NH3-NH3 

0.0000 

0.0000 

0.0026 

-0.165 

NH3-TMGa 

0.0946 

0.0945 

-0.0006 

-0.857 

NH3-TMGa:NH3 

0.0936 

0.0936 

0.0003 

-0.845 

NH3-(DMGa:NH2)3 

0.0643 

0.0643 

0.0003 

-0.646 

NH3-N2 

0.2308 

0.2314 

0.0027 

-1.820 

1.386x10'^ 

TMGa-  NH3 

0.0000 

0.0000 

0.0000 

TMGa-TMGa:NH3 

0.0271 

0.0270 

-0.0005 

0.759 

TMGa-(DMGa:NH2)3 

0.0172 

0.0172 

-0.0005 

1.116 

TMGa-N2 

0.0802 

0.0803 

0.0003 

-0.565 

TMGa:NH3-TMGa:NH3 

0.0000 

0.0000 

0.0000 

1.359x10'^ 

TMGa:NH3-(DMGa:NH2)  3 

0.0168 

0.0167 

-0.0024 

1.169 

TMGa:NH3-N2 

0.0800 

0.0800 

0.0000 

-0.567 

(DMGa:NH2)  3-(DMGa:NH2)  3 

0.0000 

0.0000 

0.0000 

4.90x10’^ 

(DMGa:NH2)3-N2 

0.0544 

0.0543 

0.0000 

-0.348 

N2-N2 

S 

0.0000 

0.0000 

0.0025 

0.151 

0.0000 

Table  6-3  Species  weight  fraction 


species 

H2 

NH3 

TMGa  TMGa;NH3 

(DMGa:NH2)  3 

N2 

COi 

8.9x10’’ 

0.376 

5.07x10’^  5.83x10’^ 

1.18x10’^ 

0.618 
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A FORTRAN  code  subroutine  Multidiffusion  based  on  these  equations  was 
developed  and  was  attached  as  Appendix,  for  future  reference.  Table  6-1  lists  Lennard- 
Jones  parameters  for  the  relevent  species.  Table  6-2  lists  the  calculated  diffusion  and 
thermal  diffusion  coefficients  of  relevent  species.  Table  6-3  lists  the  species 
concentrations  used  in  the  calculation  for  Table  6-2,  which  were  used  in  Raman 
experiment. 


Combination  Rules 

Combination  rules  for  ei2  and  oy^are  as  follow  [5,73], 


1 / 


and 
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^12  ~ ^1^2  • 6-33 

However,  when  more  physical  properties  are  available,  a more  elaborate 
formulation  may  be  used  to  achieve  up  to  10%  improvement  for  su  [80]. 

a,  =o-,[l-(C;/2.2)^  6-34 

where 


C;  6-35 

where  C^^^  is  the  dispersion  coefficient  that  describes  the  long  range  attraction.  If  C’ 

happens  to  be  greater  than  2.2,  then  the  effective  rigid-core  diameter  a\  is  taken  to  be 
zero. 


«i2  +a^) 
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then 
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C^12-«12  =^[(<^1  -ai)  + (CT2  -«2)]|l  + ^ 


ln(cTj2  -0,2)  — ln£' 
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where 


ln£:  = 


- ln(£, £2) + 3 ln(cr,  - a,  )(ct2  - «2 ) “ 


^1  -^1 

(cr,  -a,)  + (cT2  -^2) 
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and 


«]«2 

^(6) 

*^12 


a, 
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(6) 


+ - 


^2 


6-39 


where  «/  and  «2  are  the  dipole  polarizabilities,  therefore, 


^12  ~ (^1^2 ) 


(cr,  -ai)^(cr2  -^2)^ 

(cTj2  — nj2 ) 


^(6) 

'-'12 
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The  collision  integrals  were  calculated  using  the  formulation  for  a molecular  gas 
from  Reference  [80]  except  for  N2-H2,  where  the  more  elaborate  results  from  Reference 
[83]  were  used. 


CHAPTER  7 

GaN  MOCVD  SCALE  ANALYSIS  AND 
OPERATING  WINDOW  CALCULATION 

Introduction 

GaN  MOCVD  is  usually  operated  at  very  high  Peclet  number  Pe=w  JL/D,  or 
equivalently,  high  Reynolds  number  Re=wJD/v  condition  [36,37].  High  Pe  number  or  Re 
number  ensures  convective  species  transport  is  in  dominant.  As  discovered  by  in  situ 
Raman  spectroscopic  study  of  TMGa-NHa-Ni  system  and  \vill  be  demonstrated  in  the 
next  section,  outside  the  boundary  layer  at  high  Pr  or  Re  number  region,  the 
homogeneous  reaction  R-1 

TMGa+NHa  o TMGaiNHa  R-1 

should  reach  equilibrium  if  conditions  such  as  high  V/III  ratio,  high  flow  rate  and  low 
temperature  are  set.  At  these  conditions,  outside  the  boundary  layer  TMGa  concentration 
should  be  maintained  sufficiently  low  to  choke  the  oligomer  formation  reaction. 

TMGa:NH3  + TMGa+NHs  ^ (DMGa:NH2)n  R-2 

No  matter  how  high  Pe  or  Re  numbers  is,  a non-slip  boundary  condition  will  guarantee  a 
velocity  on  the  surface  equal  zero.  Therefore  a thin  boundary  layer  will  develop,  and 
inside  this  boundary  layer  it  is  always  in  diffusion-limited  regime.  Above  the  substrate 
inside  the  boundary  layer,  the  temperature  increases  rapidly,  the  ammonia  concentration 
drops  because  of  the  diffusion  limit  as  well  as  homogeneous  decomposition.  Thus,  the 
reaction  equilibrium  of  R1  will  then  breakdown. 
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Thermal  and  Concentration  Boundary  Laver  Thickness 
The  thermal  and  concentration  boundary  layer  thicknesses  in  an  impinging  flow 
are  given  as  follow  [84]; 
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dv. 


where  a is  the  strain-rate  parameter  defined  as  a = ^ , Vr=6  is  the  radial  velocity  just 

dr 

outside  the  boundary  layer,  vis  the  kinematic  viscosity,  crthe  thermal  diffusivity,  Dy  is 
the  species  diffusivity,  Sr  = v/Dy  is  the  Schmidt  number  and  Pr  = C^/u!k  is  the  Prandtl 
number,  Cp  is  the  heat  capacity  m is  the  viscosity  and  k is  the  thermal  conductivity,  a 
was  found  to  have  the  following  relation. 


a »C- 
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where  C is  a constant,  for  a stagnation  flow  over  a flat  disk  C ~ 4/;r  [85],  Ds  is  the 
diameter  of  the  susceptor  and  Uoo  is  the  velocity  at  far  away  from  the  boundary  layer.  The 
gas  phase  kinematical  viscosity  is  [86] 
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where  ref  refers  to  the  property  value  evaluated  at  a certain  reference  position. 

At  300K,  760  Torr,  the  kinematical  viscosity  of  nitrogen  and  hydrogen  are 


^-5  r 2/ 


=1.563x10  [mVs]  and  =1.095x10  ^ [mVs].  The  Prandtl  number,  Pr.  of 


N H 

nitrogen  and  hydrogen  are  ^0-713  and.  Pr^ijocA:  = 0.706 
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Using  N2  as  the  carrier  gas,  Too  =300  K,  the  thermal  boundary  layer  thickness  is 
=0.1647/ 7-9 


The  Richardson  number.  Ri,  the  relative  importance  of  buoyant  forces  to  the 


forced  convection  is  defined  as  follow, 

Ri  = {5p!  p)gUul  y_jQ 

Here  the  characteristic  length  L may  be  chosen  as  the  thermal  boundary  layer  thickness 
5x  and  Ap/p  «1 . To  overcome  the  buoyancy  effects,  Richardson  number  needs  to  be 
maintained  at  less  than  1.0.  Figure  7-1  shows  for  N2  at  Too=300K  and  760  Torr, 
Dj=0.0254  [m],  minimum  a=10  [s"'],  equivalent  to  C/oo«0.2  [m/s]  results  in  the  required 
thermal  boundary  layer  thickness  <^«0.001  [m]. 
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For  conditions  more  similar  to  a commercial  reactor,  H2  as  carrier  gas  at  inlet 
temperature  Too=300K,  substrate  temperature  Ts=1500K  and  total  pressure  P=100 
Torr,  the  expected  density  change  should  be  Ap/p  «4.  For  a reactor  with  a susceptor 
diameter  at  £>^=0.25  m,  the  minimum  required  strain  parameter  calculated  from  Equations 
7-9  and  7-10  will  be  a=5.24  [s''],  which  is  equivalent  to  the  flow  velocity  t/c»«l-0  [m/s] 
or  Rcmin  310.  The  thermal  boundary  layer  thickness  is  <^«0.019 
[m]. 


Strain  rate  [s’*] 
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Figure  7-1.  Impinging  flow  operating  window  for  Raman  experiment  {5i) 

GaN  Growth  Rate  Estimation 
No  Oligomer  Formation  Assumption 

If  the  total  pressure  is  being  set  at  a low  value,  then  the  high  order  collision 
reaction  R-2  leading  to  oligomer  formation  would  be  contained.  The  adduct  TMGa:NH3 
will  simply  decompose  and  release  TMGa,  to  under  go  homogenous  decomposition  and 
produce  MMGa  as  possible  Ga-precursor  to  arrive  to  the  growth  surface. 
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Since  the  gas  velocity  is  zero  at  the  surface,  the  molar  flux  of  the  Ga  precursor  to 
the  substrate  joa  may  be  expressed  as  [5] 


i = CnDr 
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where  Co=p/RTo  is  the  molar  gas  density  at  the  substrate  surface,  Doa  is  the  diffusion 
coefficient  for  the  Ga  precursor,  and  Xoa  is  the  Ga  precursor  mole  fraction,  and  z is  the 
distance  from  the  substrate. 


At  steady  state,  the  flux  of  Ga  precursor  to  the  substrate  will  be  equal  to  the  consumption 
rate  due  to  GaN  growth,  therefore, 

XGa  = Ror,  7-12 

Where  Ror  is  the  consumption  rate  of  Ga-precursor  due  to  GaN  growth,  given  as 


R-Gr  ~ ^Ga 


[Ga]oUc, 
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where  aoa  is  the  sticking  coefficient  of  Ga-precursor,  7/^^  is  the  mean  thermal  speed  of  a 


Ga-precursor  given  as  [5] 
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as  the  Ga-precursor  concentration  at  the  surface  of  the  substrate  is  given  as 
[Ga]o  = . 


7-16 
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At  low  Richardson  number  or  high  Pe  number  or  Re  number,  the  species  concentration 
profile  inside  the  boundary  layer  should  be  close  to  linear  [5],  therefore 
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where  , Id  can  be  view  as  a characteristic  diffusion  length  scale,  is  the  Ga-precursor 


mole  fraction  at  some  specified  reference  position,  defined  here  at  a point  outside  the 
boundary  layer.  Then  the  flux  of  the  Ga  precursor  to  the  substrate  may  be  rewritten  as 
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In  the  diffusion-limited  regime,  inside  the  boundary  layer. 
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When  substrate  temperature  To  equals  1600  K,  the  average  temperature  of 
boundary  layer  Tj  will  be  about  1400  K.  The  calculated  growth  rate,  i?or,  equals  to  1.4 
[pm/hr]  compared  with  the  experimental  result  at  1 .3  [pm/hr]  [87]  is  good  agreement. 
Including  Homogeneous  Reactions 

To  understand  homogeneous  reactions,  the  Damkohler  number  (Da),  is  introduced 
which  describes  the  homogeneous  reaction  relative  to  diffusion.  If  Daoiigomer  , the 
oligomer  formation  will  be  favored. 
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Assuming  Z=^c,  recalling  Equation  7-8,  =1.32  (Sc)  ° and 

\ ^ J \ ^ ) V'-efJ 


substitute  roo=300K,  results  in: 
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where  a again  is  the  strain  rate  with  the  unit  s'\  k\  is  a first  order  reaction  rate  constant 


with  the  same  unit. 


1)  Adduct  Formation  Reaction 


The  reaction  rate  for  adduct  formation  is  as  follow 
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where  kn  is  the  rate  constant  for  a bimolecular  collision  reaction  and  given  as 
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where  oxb  is  the  collision  effective  diameter,  /^ab  is  the  reduced  mass  and  ks  is  the 
Bolzmann  constant.  Since  the  V/III  is  greater  than  1000  and  the  ammonia  homogeneous 
decomposition  rate  is  low,  2Tnh3  is  relatively  stable  and  may  be  treated  as  a constant  to 
give  the  assumption  of  pseudo  first  order  reaction  kinetics. 
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where 
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Here  A:,  is  a pseudo  first  order  reaction  rate  constant  depending  on  both  temperature  and 


pressure. 

k_,  =k[/K^. 
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/ 0 0 0 \ 
~ ^TMGa  ~ ^NHj  ~ ^CH4  ) 

{U  — 12jMGa  ) ■ (^  ~ 


where  v=p/RT,  Kp  is  the  measured  value  given  as  3.3828xl0'^exp(AHR.3/RT)  [atm  '] 
AH=  18.5  [kcal/mol]  [58]  or  15.2  [kcal/mol]  [89].  Therefore, 
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where  k.\  is  a true  first  order  reaction  rate  constant,  is  a function  of  temperature  and 
independent  of  pressure.  The  calculated  values  of  k^  and  k.j  at  different  temperatures  are 
listed  in  Table  7-1 . Therefore,  Damkohler  number  for  adduct  formation  may  be 
represented  as: 
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and  Damkohler  number  for  adduct  decomposition  should  be 


Da.  = 


0.574aP/P. 
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When  temperature  is  higher  than  400  K,  both  k,  and  k.i  are  on  the  order  of  10 
[s'*]  as  listed  in  Table  7-1,  while  a is  on  the  order  of  10,  and  thus  local  reaction 
equilibrium  may  be  assumed. 

2)  TMGa  homogeneous  decomposition 

The  reaction  rate  for  TMGa  decomposition  is  as  follows: 
r p\ 
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where 
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3)  Oligomer  formation  reaction 
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The  reaction  rate  for  oligomer  formation  is  as  follows: 
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where  ^iMOa  is  not  constant  is  only  a pseudo  first  order  reaction  rate  constant,  A:ni  is  the 
rate  constant  for  a trimolecular  collision  reaction  given  as  [88] 


Table  7-1  shows  that  at  temperature  greater  than  700  K,  the  calculated  TMGa 
homogeneous  decomposition  reaction  rate  constant  kj  becomes  much  larger  than  the 
oligomer  formation  rate  constant,  k, , which  suggests  that  at  temperature  greater  than  700 
K,  TMGa  homogeneous  decomposition  reaction  becomes  dominant  and  the  oligomer 
formation  should  be  much  less  important.  The  TMGa  and  oligomer  concentration 
profiles  measured  by  in  situ  Raman  spectroscopic  are  shown  in  Figures  4-3  and  4-8  agree 
this  conclusion.  Therefore,  the  oligomer  formation  reaction  will  primarily  occur  at  the 
gas  phase  temperature  range  Too  < T <700K. 

Boundary  Laver  Theory  by  Karman-Pohlhausen  Method 
From  Equations  7-3  and  7-4, 
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Clearly  that  concentration  boundary  layer  is  thinner  than  thermal  boundary  layer.  As  a 
consequence,  there  will  be  a small  region  where  the  heat  transfer  is  in  the  conduction 
dominant  region  while  mass  transfer  is  still  transition  from  convection  dominated  to 
diffusion  dominated.  In  this  region,  the  gas  phase  temperature  will  increase  because  it  is 
inside  the  thermal  boundary  layer,  and  the  adduct  TMGa:NH3  will  decompose  because  of 
the  increasing  temperature.  Its  concentration,  however,  will  not  be  limited  by  diffusion 
because  it  is  still  in  the  convective  mass  transport  region. 

To  obtain  detail  information  about  this  region,  a detail  temperature  profile  inside 
the  thermal  boundary  layer  is  needed  [90].  Assuming  dimensionless  temperature©  takes 
the  form 

0 = a + 7.48) 
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where  0 = ^ , <^7-  = ^ l^e  dimensionless  thermal  boundary  layer 

Too  '^s 

thickness,  related  to  the  dimensional  boundary  layer  thickness  by  = 5j  I L . 
From  the  boundary  condition. 
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evaluating  Equation  7-48  using  Equations  7-49  and  7-50  results  in  the  following  relation. 


7-51 
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heretJj.  has  no  explicit  solution  and  can  only  be  evaluated  numerically  [90]. 

1)  Raman  Experiments 

At  the  conditions  used  in  in  situ  Raman  Spectroscopic  experiments,  the  thermal 
boundary  layer  thickness  and  concentration  boundary  layer  thickness  were  found  to  be 


=11.1  [mm]  and = 8.59  [mm].  The  ratio 


ySj  j 


= 0.774 , therefore  the  gas  phase 


temperature  at  the  edge  of  the  concentration  boundary  layer  may  be  found  by  Equation  7- 
51,  0=0.929  and  T=556.8  K.  These  calculation  results  reproduce  the  in  situ  Raman 
spectroscopic  measurement  results  in  Figures  4-4  and  4-6  very  well. 

Also  listed  in  Table  7-1  is  the  calculated  equilibrium  conversion  Xe,  the  extent  of 
reaction,  for  the  TMGaiNHs  adduct,  where  Xe=l/(1+A^e)  at  the  gas  phase  temperature  from 
T=Too  to  600K.  In  this  region  the  calculated  Xe  profile  listed  in  Table  7-1,  agrees  very 
well  with  in  situ  Raman  measurements,  indicating  that  the  concentration  dependeds  on 
temperature  only. 

2)  Commercial  Reactor  Condition 

For  the  case  of  H2  as  carrier  gas  as  discussed  in  the  first  section,  with  a inlet 
temperature  T<x)=300K,  substrate  temperature  Ts=1500K,  and  total  pressure  P=100  Torr, 
the  minimum  required  Reynolds  number,  Rcmin , was  found  to  be  3 1 0 in  order  to 
guarantee  the  Richardson  number  Ri  < 1 . Use  this  Reynolds  number  value,  the  calculated 
hydrodynamic  (^,  concentration  {5c)  and  thermal  boundary  layer  thickness  are 
0.0197  m,  0.0221  m and  0.0307  m.  Gas  phase  temperature  on  the  edge  of  hydrodynamic 

and  concentration  boundary  layer  are  Tg  ^ =504  K and  Tg^  =432  K assuming  inlet  is  at 
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least  0.0307  m away  from  the  susceptor.  From  Table  7-1,  oligomer  formation  below 
these  temperatures  inside  this  transition  layer  may  be  kept  at  a minimum. 

The  Damkohler  number  for  oligomer  formation  may  be  written  as 
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Assuming  local  equilibrium  of  Rl,  the  concentration  of  TMGa  may  be  written  as 
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where  is  the  average  mole  fraction  of  TMGa  in  the  temperature  from  Too  to  600K. 
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where  is  the  reaction  equilibrium  constant  calculated  at  the  average  temperature  T 
Assuming  average  temperature  equals  500  K. 
k" 
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The  opterating  window  is  hightlited  in  Figure  7-2  where  boundary  layer  thickness  is 
slightly  large  than  2 [mm]. 
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Figure  7-2.  GaN  Operating  Window  (N2) 


Table.  7-1  various  homogeneous  reaction  rate  constants  at  atmospheric  pressure 


Temp 

k-i 

k,’ 

kt 

ke 

k," 

300 

4.8E+04 

1.9E+09 

1.5E-29 

4.0E+04 

5.8E-04 

400 

2.4E+07 

1.6E+09 

l.OE-18 

6.8E+01 

3.6E-04 

500 

l.OE+09 

1.5E+09 

3.4E-12 

1.4E-00 

9.7E-05 

600 

l.lE+10 

1.3E+09 

7.4E-08 

l.lE-01 

2.2E-06 

700 

6.7E+10 

1.2E+09 

9.2E-05 

1.8E-02 

5.6E-08 

800 

2.4E+11 

l.lE+09 

1.9E-02 

4.8E-03 

3.0E-09 

900 

6.7E+11 

l.lE+09 

1.2E+00 

1.6E-03 

3.0E-10 

1000 

1.4E+12 

l.OE+09 

3.4E+01 

7.1E-04 

4.8E-11 

1100 

2.8E+12 

l.OE+09 

5.2E+02 

3.5E-04 

l.OE-11 

1200 

4.8E+12 

9.6E+08 

5.0E+03 

1.9E-04 

2.8E-12 

To  revisit  the  example  by  recalling  Equation  7-25: 
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When  the  substrate  temperature  To=1400  K,  the  average  temperature  of  boundary  layer 
Tg  will  be  1 000  K,  transition  region  temperature  T =600K.  T calculated  Damkohler 

number  equals  0.004  and  the  calculated  rate  of  growth  Ror  = 1.408  x(l-Da)  « 1.4  qm/hr. 

Additional  Comments 

Growth  Rate 

It  has  been  widely  reported  that  growth  rate  of  GaN  is  always  in  the  rage  of  1 -4 
[pm/hr],  regardless  of  the  conditions.  Inside  thermal  boundary  layer  outside  the 
concentration  boundary  where  adduct  and  TMGa  concentrations  should  depend  on 
temperature  only.  The  hydrodynamic  boundary  layer  thickness  from  Karman-Pohlhausen 
method  [90]  was  found  to  be 

= 1.3751  Re 7-56 
For  the  conditions  used  in  GaN  MOCVD, 

S < S(j  < Sj.  7-57 

If  the  gas  phase  temperature  at  the  edge  of  the  hydrodynamic  boundary  layer  is  about  600 
K,  the  mole  fraction  of  the  adduct  according  to  Table  7-1  will  be  always  lower  than  ten 
percent  at  the  edge  of  the  concentration  boundary  layer.  Therefore  the  Damkohler 
number  of  the  oligomer  formation  inside  the  hydrodynamic  boundary  should  always  be 
small.  Most  oligomer  was  formed  outside  the  hydrodynamic  boundary  layer  and  thus 
carried  away  by  the  main  flow  resulting  low  Ga  yield  but  steady  growth  rate. 

Pressure  dependence 
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Outside  the  boundary  layer  in  the  transitional  region,  the  temperature  should  be 
lower  than  500  K,  the  p times  accoding  to  the  calculation  results  listed  in  Table  7-1  is 
always  much  greater  than  1 . Therefore,  Equation  7-54  becomes: 
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Da  is  independent  of  pressure. 

Inside  the  boundary  layer,  the  temperature  should  be  higher  than  500K,  p times  Kp 
becomes  smaller  than  1.  Therefore,  Equation  7-54  becomes: 


Da-^^^Na^KG^ 


a 


ABC 


' SkgT ' 

U.5 

( \ 

P 

RT, 

Ga 


7-58b 


Da  is  a function  of  p^.  Pressure  has  very  strong  influence  on  oligomer  formation. 
Carrier  Gas  Effect 


Because 


Ni 

« 0.657  and 

^ c 

\^T  y 

y 

0.434,  using  H2  as  the  main  carrier  gas 


will  surely  inducing  broader  transition  region  resulting  in  a higher  level  of  oligomer 
formation. 

Two-Flow  Configuration 

In  the  two-flow  configuration,  H2  is  chosen  as  the  main  carrier  gas  that  isdirectly 
introduced  into  the  thermal  boundary  layer,  secondary  flow  of  N2  gas  flows  outside  the 
bounday.  Recalling  Equation  7-25, 
k. 


Da  = 


DabIL' 


7-25 


Outside  the  boundary,  the  N2  gas  will  reduces  the  transition  region  and  therefore,  reduce 
the  oligomer  formation.  Inside  the  boundary,  under  the  assumption  of  pseudo  first  order 
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reaction,  the  rate  constant  ki  remains  the  same.  Thus  using  a main  carrier  gas  with  large 
diffusivity  such  as  H2  will  further  decrease  the  oligomer  formation. 


CHAPTER  8 

NUMERICAL  SIMULATION  OF  TMGa-NH3-N2  FLOW-REACTION  SYSTEM 


Introduction 


A mathematical  model  of  MOCVD  of  GaN  using  TMGa  and  NH3  that  is  built 
upon  the  fundamental  transport  phenomena  and  assumed  reaction  mechanisms  have  been 
used  by  several  research  groups  [9,1 1,63-66]  to  provide  accurate  simulation  of  fluid 
pattern,  temperature  and  species  concentration  profiles  for  reactor  design  and  process 
optimization.  No  in  situ  experimental  measurements  on  the  gas  phase  temperature  and 
species  concentration  have  been  available  for  comparison  to  validate  the  model.  Based 
on  the  formulation  for  mass  and  thermal  diffusion  coefficients  that  are  discussed  in 
Chapter  6 and  the  reaction  mechanism  proposed  in  Chapter  4,  a new  species 
mathematical  model  is  used  to  simulate  the  fluid  dynamics  and  chemical  reactions  of  the 
TMGa-NH3-N2  system.  The  simulation  results  will  be  used  to  compare  with  the  in  situ 
Raman  spectroscopy  experiment  results. 

The  number  of  molecules  of  species  j can  be  calculated  from  its  relative  Raman 
cross  section  using  Equation  3-15: 


Ij  (da/dn)j  Nj 


3-15 


(dcjIdQ),/ 


Therefore, 


Ij  f (daldQ) 


8-1 


For  a binary  system,  the  mole  fraction  of  species  j is 
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“jh . 

■ I j (Scr/aQ)^  ’ 

Jn,  {dcrldOL)^^ 

and  for  multicomponent  system,  the  mole  fraction  of  species  j is 


Ij  ( {dc7ldO)j 

In,  [(ao-/5Q)^^ 


J 

1+AJ 

' (a<T/aQ),  ' 

-1 

, I2 

' (dcj/dn)^  ^ 

-1 

' (da/dQ)^  ' 

V 

(acr/SQ)^^  ^ 

jacr/aQ)^^  ^ 
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Both  the  absolute  and  relative  Raman  cross  sections  of  the  adduct  TMGaiNHa  and 
oligomer  (DMGa:NH2)x,  are  not  know. 

From  Equations  8-2  and  8-3,  the  species  relative  Raman  scattering  cross  sections 
can  not  be  separated  from  the  concentration  terms,  therefore,  in  a flow  system,  direct 
identification  of  the  relative  Raman  cross  section  of  species  j from  the  measured  Raman 
intensity  ratio  is  mathematically  forbidden.  Any  such  attempt  will  only  result  in  a 
arbitrary  number  except  in  a binary  system  at  room  temperature,  in  which  physical 
properties  such  as  diffusivity  and  viscosity  are  known,  and  the  species  relative  Raman 


scattering  cross  section 


{dcrldn)j 

{dcrl 


is  the  only  unknown  constant  parameter,  will  the 


identification  result  in  a unique  solution. 

Because  key  species  physical  properties  such  as  Raman  cross  section  and  vapor 
pressure  of  the  inlet  concentration  values  of  some  key  components  were  not  available, 
simulation  was  conducted  with  assumed  values. 
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Experimental 

The  gas  phase  reaction  system  equipped  with  in  situ  Raman  spectroscopy  as 
described  in  Chapter  2 was  used  in  this  study  (Figures  2-5,  2-6  and  4-1).  High  purity  N2 
(99.999%)  at  a flow  velocity  of  3.0  cm/s  was  introduced  into  the  center  tube,  with  50% 
weight  fraction  of  high  purity  “blue”  NH3  (99.999%)  in  N2  carrier  gas.  In  a separate  line 
high  purity  N2  (99.999%)  was  bubble  through  a TMGa  bubbler  to  transport  TMG.  This 
flow  was  combined  with  the  NH3/N2  stream  to  yield  a NHs/TMGa  (V/III)  mole  ratio 
equal  to  500.  The  center  flow  tube  was  packed  with  3mm  diameter  quartz  beads.  A 
second  flow  of  high  purity  N2  (99.999%)  was  delivered  to  the  oter  regin  at  the  same 
average  velocity.  The  flow  velocity  was  calculated  at  standard  conditions  of  296.15  K 
and  one  atmospheric  pressure.  The  reactor  aspect  ratio  was  adjusted  to  1 .0.  The  heater 
temperature  was  set  at  1373.15  K measured  by  a thermocouple  and  the  reactor  pressure 
was  kept  at  one  atmosphere.  Raman  spectra  were  taken  along  the  center  line  of  the 
reactor.  The  species  identification  was  discussed  in  Chapter  4. 

N2  rotational  part  of  the  Raman  spectra  was  used  for  gas  phase  temperature 
measurement,  while  N2  at  2331  cm'*,  NH3  at  3337  cm"',  TMGa  at  526  cm'',  TMGa:  NH3 
at  528  cm"'  and  the  oligomer  524  cm''  vibrational  Raman  transition  peaks  were  used  for 
species  concentration  profile  measurements.  Quantitative  interpretation  of  TMGa, 
TMGa:  NH3  and  oligomer  require  a sophisticated  deconvolution  procedure.  Details  of 
the  data  reduction  procedure  are  described  in  Chapters  3 and  4.  Gas  phase  temperature 
and  species  concentration  profiles  along  the  center  line  of  the  reactor  are  shown  in 
Figures  4-4,  4-5  and  4-3. 
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Gas  Phase  TMGa-NHvN?  Simulation 

Governing  Equations 

The  model  is  a axisymmetric  one  that  solved  the  steady  state  momentum, 
energy  and  mass  conservation  equations  in  cylindrical  coordinates  [5]. 
r-direction 


dVf. 

P(vr 

dr 
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dz 
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dz 


Energy  conservation  equation: 

\ d dT^  5^1  . ^ 

ir — ) + — T 5-9 

r dr  dr 

Species  equations  of  change 

- V/7, -V- V- j, +r,  = 0 6-10 

and 

j =-p,bJV\nT — -«,Ve)  6-12 

Gas  Phase  Chemical  Reaction  Model 

Gas  phase  chemical  reactions  are  discussed  in  Chapter  4,  reactions  G1-G4  are 
proposed  from  the  results  of  in  situ  Raman  experiments  discussed  in  that  chapter. 


pc  ivr—  + vz—)  = k 
P dr  dz 
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The  reaction  rate  constants  were  discussed  in  Chapter  6 and  listed  in  Table  8-1. 

Table  8-1.  Estimated  gas  phase  reaction  rate  constants 
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Galerkin  Finite  Element  and  Consistent  Penalty  Method 

A finite  element  computational  method  for  two-dimensional  laminar 
incompressible  flows  is  described  below.  The  method  is  based  on  the  consistent  penalty 
finite  element  method  [91]. 

The  Navier-Stokes  equations  are  given  as  [5]: 
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in  Q 


8-4 


in  Q 


8-5 


dX  : 


where  Q is  the  open  bounded  domain,  the  subscripts  i and  j denote  the  coordinate 
directions,  p is  the  density,  u,  is  the  velocity  component  in  the  i-th  coordinate  direction,  p 
is  the  pressure,  p is  the  molecular  viscosity  of  the  fluid,  b,  is  the  body  force  in  the  i-th 
coordinate  direction,  and  5,y  is  the  Kronecker  delta. 

In  the  penalty  method,  the  conservation  of  mass  equation  (the  continuity  equation) 
is  expressed  as: 


In  the  weighted  residuals  method  [86],  the  test  function  for  the  momentum 
equation  and  the  continuity  equation  are  denoted  as  Wu(x)  and  Wp(x).  Therefore,  the 
weak  forms  of  the  Navier-Stokes  equations  are  as, 


5u,  _ 1 

dx , 
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to  eliminate  the  pressure  term,  where  A,  is  the  penalty  number. 
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Integrated  by  parts,  Equation  8-7  becomes 
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Using  finite  element  discretization  notation,  integrated  over  entire  flow  domain  Q, 
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dx  = 0 


8-11 


In  the  Galerkin  finite  element  method,  the  test  functions  are  selected  from  the 


same  space  of  interpolating  polynomials  as  the  trial  functions.  Therefore, 


where 


C = [ {pVCnI  9);^dx 

•“e-  " - aXj. 
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8-19 

Uj  is  a column  vector  of  nodal  values  of  the  velocity  component  Ui,  p is  a column  vector  of 
nodal  pressure  f is  a column  vector  of  interpolating  polynomials  for  velocity,  s is  a 
column  vector  of  interpolating  polynomials  for  pressure,  {b.c.}  is  a column  vector  of  the 
flux  boundary  condition,  and  the  subscripts  i and  j denote  the  spatial  dimensions.  The 
second  order  stress  tensor  term  is  integrated  by  parts  using  the  Green-Gauss  theorem;  the 
continuous  flow  variables  are  interpolated  using  the  nodal  values  of  these  variables  and 
the  appropriate  interpolation  polynomials;  the  weak  form  Navier-Stokes  equations  are 


Figure  8-1 . Two-dimensional  flow  (square  9 nodes)  and  pressure  (triangle  3 nodes) 
elements 

integrated  over  each  element  to  obtain  the  discrete.  The  integrations  in  Equations  8-14 
through  Equation  8-19  were  evaluated  using  the  Gauss  numerical  quadrature  method  with 
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five  Gauss  points  in  each  coordinate  direction.  Figure  8-1  shows  the  two-dimensional 
flow  (square)  and  pressure  (triangle)  elements. 

Results  and  Discussion 

The  simulated  gas  phase  temperature  profile  along  the  center  line  of  the  reactor  is 
shown  in  Figure  8-1.  The  profile  seams  to  fit  the  Raman  experimental  results  of  Figure 
4-4  very  well. 


Figure.  8-1  Gas  Phase  Temperature  profile 

The  simulated  NH3  concentration  profile  at  the  same  position  is  shown  in  Figure 
8-2.  It  also  fits  the  in  situ  Raman  experimental  results  of  Figure  4-5  very  well.  This 
comes  with  no  surprise  because  NH3  and  carrier  gas  N2  composed  more  than  99%  of  the 
gas  phase,  and  their  properties  are  well  known.  Besides,  NH3  decomposes  very  little  both 
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homogeneously  and  heterogeneously  and  their  reaction  rate  constants  were  well 
identified  in  Chapter  5. 


Figure  8-2.  Simulated  NH3  weight  fraction  profile. 

The  simulated  TMGa,  adduct  and  oligomer  concentration  profiles  are  shown  in 
Figures  8-3,  8-4  and  8-5.  Comparing  with  the  in  situ  Raman  experimental  results  shown 
in  Figure  4-8,  the  simulation  produces  excellent  resemblance  of  the  experimental  TMGa, 
adduct  and  oligomer  concentration  profiles.  This  indicates  that  the  gas  phase  reaction 
mechanisms  and  their  reaction  rate  constants  proposed  in  Chapters  4 and  7 provide  a 
realistic  approximation.  The  simulation  successfully  predicted  the  rapid  disappearence  of 
the  TMGa  at  several  mm  from  the  reactor  inlet.  TMGa  reapperes  and  increase 
dramatically  close  to  the  thermal  boundary  layer  because  of  the  disassociation  of  the 
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adduct.  Inside  the  thermal  boundary  layer,  with  the  steep  gas  phase  temperature  gradient, 
TMGa  decomposes.  Its  concentration  after  peaked  at  about  700  K,  starts  to  descent 
again. 


0 5 10  15  20  25  30  35 


Distance  from  the  Suscepter  [mm] 

Figure  8-3.  Simulated  TMGa  weight  fraction  profile 

Simulation  also  predicts  the  fast  descent  of  the  adduct  concentration  at  the  region 
close  to  the  inlet,  where  high  TMGa  concentration  promotes  fast  oligomer  formation. 
After  the  early  descent,  the  adduct  concentration  remain  relatively  constant  as  does  the 
gas  phase  temperature  as  shown  in  Figure  8-1 . Here  the  inlet  TMGa  is  mostly  depleted, 
since  gas  phase  temperature  is  not  sufficiently  high  to  prompt  the  adduct  dissociation; 
therefore,  an  equilibrium  state  has  been  reached.  Also,  thermal  diffusion  plays  an 
important  role  in  keeping  the  concentration  of  the  large  molecule  such  as  the  adduct  and 
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the  oligomer  in  the  relatively  low  temperature  region. 


Figure  8-4.  Simulated  TMGa:NH3  weight  fraction  profile 

Upon  entering  the  thermal  boundary  layer,  the  gas  phase  temperature  begins  to 
rise  and  the  adduct  begins  to  deassociate  form  TMGa  and  NH3.  The  TMGa  then  reacts 
with  the  remaining  adduct  to  form  Oligomer.  This  causes  the  concentration  of  the  adduct 
to  decrease  dramatically  in  the  vicinity  and  inside  the  thermal  boundary  layer.  On  the 
other  hand,  the  oligomer  concentration  as  shown  in  Figure  8-5  experiences  two  drastic 
increase  stages;  one  at  the  region  close  to  the  inlet  where  TMGa  concentration  still  high, 
and  another  in  the  vicinity  of  the  thermal  boimdary  layer  where  the  TMGa  concentration 
surges  when  gas  phase  temperature  begins  to  rise.  The  simulation  predicts  the  fast 
decrease  of  the  oligomer  concentration  at  near  the  suscepter  region,  which  is  the  only 
discrepancy  that  can  be  found  between  the  simulation  results  and  Raman  experimental 
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results. 


Figure  8-5.  Simulated  Olegomer  weight  fraetion  profile. 

There  are  everal  possible  causes  that  may  contribute  to  this  discrepancy.  First,  the 
adduct  and  the  oligomer  have  very  low  vapor  pressures.  They  may  physically  condense 
[94]  in  the  low  temperature  to  alter  their  gas  phase  concentration  profiles.  This 
phenominum  was  not  included  in  this  simulation. 

Second,  the  simulation  was  severely  restricted  by  the  computational  copacity. 
Currently,  the  computer  can  only  support  the  simulation  of  two  velocity  components,  one 
temperature  and  six  species  concentration,  for  a total  of  nine  partial  differencial 
equations.  Many  small  molecules  such  as  dimethylgallium  (DMGa),  monomethylgallium 
(MMGa),  methane  (CH4)  and  NH2  cannot  be  included  in  this  simulation.  Their 
concentration  will  be  made  up  by  carrier  gas,  N2.  Since  NFls  homogeneously 
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decomposes  to  a small  extent  and  a high  V/III  ratio  was  used  in  the  experiment,  NH2 
concentration  may  still  be  significant  when  comparared  with  that  of  the  oligomerand.  At 
the  same  time,  TMGa  completely  decomposes  into  smaller  molecules.  Therefore,  DMGa 
and  MMGa  concentration  could  be  high  when  compared  with  that  of  the  oligomer. 
Therefore,  N2  concentration  may  be  over  estimated  in  a relatively  large  extent. 


Temperatue  [k] 
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Figure  8-6.  Simulated  flow  steam  line,  and  spatial  distribution  of  gas  phase  temperature 
and  TMGa  concentration  of  the  entire  flow  region 

Inside  the  thermal  boundary,  the  gas  phase  temperature  gradient  is  very  high. 

This  temperature  gradient  will  drive  the  heavy  molecules  away  from  the  hot  substrate 

surface,  at  the  same  time  pull  small  molecules  into  the  thermal  boundary  layer  close  to 
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the  substrate  surface.  The  over  estimation  of  light  molecule  concentration  will  centainly 
cause  the  simulation  over  estimate  the  thermal  diffusion  effect  of  the  heaviest  molecule 
the  oligomer. 

Third,  the  region  that  shows  the  discrepancy  is  close  to  the  suscepter,  where  the 
suscepter  blocked  a significant  portion  of  the  Raman  signal  and  thus  increases  the  error. 
In  addition  the  high  gas  phase  temperature  in  this  region  generates  higher  nosie  level 
which  further  uncertainty  in  the  experimental  results. 
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Figure  8-7.  Simulated  TMGa,  adduct  and  oligomer  concentration  of  the  entire  flow 
region 
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Of  cause,  possible  unidentified  species  may  be  formed  in  the  high  temperature 
region  and  their  Raman  signal  eould  overlap  with  those  of  the  speeies  to  produce  error  in 
the  eoncentration  ealculations. 

The  simulated  flow  stream  line,  and  spatial  distribution  of  gas  phase  temperature 
and  NH3  concentration  of  the  entire  flow  region  are  shown  in  Figure  8-6.  The  TMGa, 
adduct  and  oligomer  concentrations  of  the  entire  flow  region  are  given  in  Figure  8-7. 

Conelusions 

Numerical  simulation  of  TMGa:NH3-NH3-N2  gas  phase  reaction  system  was 
successfully  carried  out.  Simulated  gas  phase  temperature  and  NH3  eoncentration 
profiles  fit  the  in  situ  Raman  experimental  results  very  well.  The  simulated  eoncentration 
profiles  of  other  species  also  produeed  realistic  resemblance  of  the  in  situ  Raman 
experimental  results.  Laek  of  species  vapor  pressure  and  Raman  cross  section  data 
prevent  the  identification  of  reaction  rate  constants.  However,  the  theoretical  estimation 


produees  good  results. 


CHAPTER  9 

CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FUTURE  WORK 

Conclusions 

Method 

Metal  Organic  Chemical  Vapor  Deposition  (MOCVD),  is  a convoluted  process  of 
several  individual  processes  such  as  heat  transfer;  species  mass  transfer,  homogeneous 
chemical  reactions,  and  heterogeneous  chemical  reactions: 

1-1 

When  knowledge  of  the  three-dimensional  spatially  resolved  gas  phase  temperature  (7), 
species  concentration  (x)  and  reaction  mechanisms  {M)  have  been  obtained,  operators  F 
and  R can  be  uniquely  determined,  the  mappings  of  convolution  and  the  deconvolution 
become  isomorphic  [15]. 

Go[F(y,r,P,x)]©[/?(r,P,x,M)]  1-12 

By  using  in  situ  Raman  spectroscopy,  the  gas  phase  temperature  (7),  species 
concentration  (x)  and  reaction  mechanisms  {M)  along  the  centerline  of  the  reactor  were 
determined.  These  experimental  data  were  then  subjected  to  reactor  inverse  modeling  for 
the  deconvolution  of  G into  F and  R by  which  “unknown”  physical  and  chemical 
properties  of  some  gas  phase  species  may  be  identified. 

Relative  Normalized  Differential  Raman  Scattering  Cross  Section  Modification 

To  accurately  determine  species  concentration  at  high  gas  phase  temperature 
using  in  situ  Raman  spectroscopy,  modification  of  the  "Relative  normalized  differential 
Raman  scattering  cross  section  “Sj"  formula  was  proposed. 
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1 - exp(-2331cw  ^hc/kT) 
1 - exp(-3337cm~^/?c/A:7’) 


3-23 


Parameter  Identification 

For  the  first  time  the  gas  phase  binary  mass  diffusion  coefficient  of  TMGa  in  N2 
was  successfully  identified.  The  gas  phase  binary  mass  diffusion  coefficient  of  TMGa 
and  N2  was  found  to  be  0.08  [cm  s'  ]. 

The  NH3  homogeneous  decomposition  rate  constant  was  successfully  identified. 
The  reaction  rate  constant  of  the  first  step  ki  of  NH3  homogeneous  decomposition 
reaction  rate  constant  was  found  to  be  1.4xl0’^exp(90.1  IkCal/RT)  [cm^  mof'  s'*]. 

TMGa  homogeneous  decomposition  rate  constant  was  also  successfully 
identified.  The  TMGa  decomposition  reaction  rate  constant  ki  was  found  to  be  3.4x1 0'^ 
exp(60.1  kCal/RT)  [s'*]. 

The  success  in  finding  a unique  solution  to  the  two  parameters  of  a binary  mass 
conservation  equation  suggests  that  using  one  single  boundary  measurement,  the 
centerline  measurement,  is  sufficient  for  uniquely  identifying  these  parameters  in  a binary 
system.  For  the  binary  system,  the  problem  is  well  posed. 

Reaction  Mechanism  Identification 

The  unique  properties  of  in  situ  Raman  spectroscopy  can  provide  3-D  spatially 
resolved  gas  phase  temperature  and  species  concentration  have  also  proven  to  be  useful 
for  understanding  gas  phase  reaction  mechanisms.  For  the  first  time  TMGa-NH3-N2 
reaction  system  was  studied  at  conditions  close  to  the  real  GaN  MOCVD  conditions.  The 
high  order  reactions  path  can  be  studied  at  relatively  high  reaction  pressure.  Gas  phase 


122 


reaction  mechanisms  were  proposed  based  on  the  species  concentration  profiles  of  the 
key  components  of  the  TMGa-NH3-N2  reaction  system, 

ki 

NH3  + M ► NH2  + H G-1 

ki  k2  k3 

TMGa ► DMGa  ► MMGa ► Ga  G-2 


TMGa+NH3  ■ ■>  TMGa:NH3 


G-3 


TMGa:NH3  + TMGa  ► DMGa:NH2:TMGa  + CH4 

DMGa:NH2:TMGa  + NH3  ► (DMGa:NH2)2  + CH4  G-4 

A multicomponent  mathematical  CVD  model  was  developed  to  numerically 
simulate  the  gas  phase  temperature  and  species  concentration  profiles  using  the  proposed 
reaction  mechanisms.  The  simulated  results  fit  the  in  situ  Raman  results  very  well. 

GaN  MOCVD  Process  Scale  Analysis 

With  the  proposed  TMGa-NH3-N2  reaction  mechanism,  scale  analysis  and 
operating  window  calculation  for  GaN  MOCVD  provided  qualitative  understanding  of 
MOCVD  of  GaN.  The  analytical  equation  for  estimating  the  GaN  growth  rate  was 
deduced. 


_ xl.7 


T 


l32RT„(pJ’{uJ 


-WTXSi  (l-Oa) 


7-25 


where  Da  is  the  Damkohler  number  for  oligomer  formation 
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Da  = 


0.514aP/P^^f 


^ l + p-Kp 


Adduct 


7-54 


An  optimized  operating  window  for  GaN  growth  was  also  suggested.  Both  agree  very 
well  with  the  operating  conditions  provided  by  commercial  reactor  manufacturers. 
Numerical  simulation  of  TMGa:NH j-NH^-N?  gas  phase  reaction  system 

Numerical  simulation  of  TMGa:NH3-NH3-N2  gas  phase  reaction  system  was 
successfully  carried  out.  Simulated  gas  phase  temperature  and  NH3  concentration 
profiles  fit  the  in  situ  Raman  experimental  results  well.  Simulated  concentration  profiles 
of  other  species  also  produced  realistic  resemblance  of  the  in  situ  Raman  experimental 
results.  Lack  of  species  vapor  pressure  and  Raman  cross  section  data  prevent  the 
identification  of  reaction  rate  constants.  However,  the  theoretical  estimation  produces 
good  results. 

Table  8-1 . Estimated  gas  phase  reaction  rate  constants 


NH3+M 


TMGa+  NH3 


TMGa 


TMGa:NH3  ^-1 


IC^  _ ^11 


pKp  KpRT 


TMGa:NH3 
+ TMGa  + NH3 


0.5 
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Recommendation  for  Future  Work 

The  combination  of  Raman  spectroscopy  and  reactor  inverse  modeling  has  proven 
to  be  effective  for  obtaining  both  physical  and  chemical  properties  of  gas  phase  species  in 
a binary  system.  For  a multicomponent  system,  it  is  not  clear  whether  the  problem  is 
well  posed. 

• Does  a single  boundary  measurement  contain  sufficient  information  for  a unique 
solution? 

• Which  solution  structure  Dirichlet-to-Neumann  mapping  or  others  can  guarantee 
the  existing  of  a unique  solution? 

• Which  solution  method  can  guarantee  locating  an  unique  solution? 
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APPENDIX 

SUBROUTINE  MULTIDIFFUSION 


############################################################### 

SUBROUTINE  MULTIDIFFUSION 

###################################################  DECLARATION 
IMPLICIT  REAL* 8 (A-H,0-Z) 
common/temp/ at(20) 

common/integrl/omijl  1(55, 20), oml2ij(55, 20), omij  13(55,20), 

omij22(55,20) 

common/spsinx/nofsps,non,np(  1 0, 1 0),mp(  10,10) 
common/species/  vmi(10),yi(10,20),xi(10,20),epkij(55),sigij(55) 
common/dij  dti/dij  1 (5  5 ,20),dij2(5 5 ,20),bdij  (5  5 ,20),dti(5 5 ,20) 
dimension  vmm(20),omigal  1 (20),omiga  1 2(20),omiga  1 3 (20),omiga22(20) 

open  (unit=l,  name="out.dat"  ,status='new') 

non=l 

do  i=l,non 

at(i)=300. 

end  do 

nofsps=6 

call  diffusion 

write  (1,*)  "dijl(kj,i),dij2(kj,i),bdij(kj,i)" 
do  i=l,non 
kj=0 
sigb=0.0 
do  k=l,nofsps 

do  j=k,nofsps 
kj=kj+l 

write  (1,*)  dijl(kj,i),dij2(kj,i),bdij(kj,i) 

end  do 

end  do 

write  (1,*)  "yi(k,i),bdij(kj2,0" 
do  k=l,nofsps 
J2=l 

kj2=np(k,j2) 

sigb=sigb+bdij(kj2,i)*yi(k,i) 

print  *,  "yi(k,i),bdij(kj2,i)",yi(k,i),bdij(kj2,i) 


n n 
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write  (1,*)  yi(k,i),bdij(kj2,i) 

end  do 

print  *,  "sig  lu-i*bdi-l=",sigb 
write  (1,*)  "sig  lu-i*bdi-l=",sigb 
write  (1,*)  "dti(k,i),lu*dt(k,i)" 
sigdt=0.0 
do  k=l,nofsps 

write  (1,*)  dti(k,i),yi(k,i)*dti(k,i) 
sigdt=sigdt+yi(k,i)*dti(k,i) 

end  do 

write  (1,*)  "sig-lu-i*dT-i=",sigdt 


C ################################################### 


C ################################################### 

c 

SUBROUTINE  MATOPR2  (OPERAT,vMATl,vMAT2,vMAT3,NORDER) 

###################################################  DECLARATION 
IMPLICIT  REAL*  8 (A-H,0-Z) 

CHARACTER* 003  OPERAT 

DIMENSION  IVEC(50),vMAT  1(10,1 0),vMAT2(  1 0, 1 0),vMAT3(  1 0,10) 

C #######################################  CHOOSE  MATRIX 

OPERATION 

write  (1,*)  "something" 

IF  (OPERAT  .EQ.  'ADD')  GO  TO  10 
IF  (OPERAT  .EQ.  'IDN')  GO  TO  20 
IF  (OPERAT  .EQ.  'MUL')  GO  TO  30 
IF  (OPERAT  .EQ.  'SUB')  GO  TO  40 
IF  (OPERAT  .EQ.  'INV')  GO  TO  50 

C ###############################################  MATRIX  ADDITION 

10  DOI=l,NORDER 
DO  J=l,NORDER 

vMAT3(I,J)  = vMATl(I,J)+vMAT2(I,J) 

END  DO 
END  DO 
RETURN 

C ###############################################  IDENTITY  MATRIX 

20  DOI=l,NORDER 
DO  J=l,NORDER 
IF  (I  .EQ.  J)  vMATl(I,J)=1.0 


end  do 
close  unit=l 
end 


C ## 

C ## 

c ## 


SUBROUTINE  MATOPR2 
FEM  VERSION  2.0 


min  huang 


## 

## 

## 
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IF  (I  .NE.  J)  vMATl(I,J)=0.0 
END  DO 
END  DO 
RETURN 

c #########################################  matrix 

MULTIPLICATION 
30  DOI=l,NORDER 
DO  J=  UNORDER 
vMAT3(I,J)  = 0.0 
DO  K=  UNORDER 

vMAT3(I,J)  = vMAT3(I,J)+vMATI(I,K)*vMAT2(K,J) 

END  DO 
END  DO 
END  DO 
write  (U*)  "mul" 

RETURN 

c matrix 

SUBTRACTION 
40  DO  1=  unorder 

DO  J=  UNORDER 

vMAT3(I,J)  = vMATI(I,J)-vMAT2(I,J) 

END  DO 
END  DO 
RETURN 

c mmmmummmmmmmmmmmmummmmm  matrix  inversion 

50  DO  I- unorder 

DO  J=  UNORDER 
vMAT3(I,J)  = vMATl(I,J) 

END  DO 
END  DO 

DO  1=  UNORDER 
IVEC(I+20)  = I 
END  DO 

DO  1=  unorder 
R1  =0.0 
M =I 

DO  J=I,NORDER 

IF  (dABS(Rl)  .LT.  dABS(vMAT3(I,J)))  THEN 
M =J 

RI  = vMAT3(U) 

END  IF 
END  DO 

IF  (I  .NE.  M)  THEN 
K = IVEC(M+20) 

IVEC(M+20)  = IVEC(I+20) 
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IVEC(I+20)  = K 
DO  J=l,NORDER 
S = vMAT3(J,I) 
vMAT3(J,I)  = vMAT3(J,M) 
vMAT3(J,M)  = S 
END  DO 
END  IF 

vMAT3(I,I)=  1.0 
DO  J=l,NORDER 
vMAT3(I,J)  = vMAT3(I,J)/Rl 
END  DO 

DO  J=l,NORDER 
IF  (I  .EQ.  J)  GO  TO  60 
R1  = vMAT3(J,I) 

IF  (dABS(RI)  .LE.  l.OE-20)  GO  TO  60 
vMAT3(J,I)  = 0.0 
DO  K=I,NORDER 

vMAT3(J,K)  = vMAT3(J,K)-Rl*vMAT3(I,K) 

END  DO 

60  END  DO 

END  DO 

DO  I=I,NORDER 
IF  (IVEC(I+20)  .EQ.  I)  GO  TO  90 
M = I 

70  M = M+1 

IF  (IVEC(M+20)  .EQ.  I)  GO  TO  80 
IF  (NORDER  .GT.  M)  GO  TO  70 

80  IVEC(M+20)  = IVEC(I+20) 

DO  J=  I, NORDER 
R1  = vMAT3(I,J) 
vMAT3(I,J)  = vMAT3(M,J) 
vMAT3(M,J)  = R1 
END  DO 
IVEC(I+20)  = I 
90  END  DO 

write  (1,*)  "inv" 

RETURN 

C ###########################################################  END 

END 

C ################################################### 

C ##  SUBROUTINE  diffusion  ## 

C ##  FEM  VERSION  2.0  ## 

C ##  min  huang  ## 

C ################################################### 

c 


o n 
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SUBROUTINE  diffusion 

declaration 

IMPLICIT  REAL* 8 (A-H,0-Z) 

common/integrl/omij  11(55 ,20),omij  12(5  5 ,20),omij  1 3 (5  5 ,20), 

& omij22(55,20) 

common/spsinx/nofsps,non,np(  1 0, 1 0),mp(  10,10) 
common/species/  vmi(10),yi(10,20),xi(10,20),epkij(55),sigij(55) 
common/dijdti/dijl(55,20),dij2(55,20),bdij(55,20),dti(55,20) 
common/temp/ at(20) 

dimension  vmm(20),omigal  I(20),omigal2(20),omigal3(20),omiga22(20) 


& ,vm(10),sigii(10),epii(10) 

c - find  i 


C 


c 

h2 

11  12 

13 

14 

15  16  1 2 

3 

4 

5 

6 

c 

nh3 

10  11 

22 

23 

24  25  26 

7 

8 

9 

c 

tmg 

14  15 

33  34  35  36 

12 

13 

c 

tmg:nh3 
16  17 

18 

44  45  46 

c 

oligomer 
19  20 

55  56 

c 

n2 

21 

66 

c 

c 

l:h2 

7:nh3 

12:tmg  16:tmg-nh3 

19:oligomer 

21:n2 

c 

c pointer 

np(l,l)=l 

np(l,2)=2 

np(l,3)=3 

np(l,4)=4 

np(l,5)=5 

np(l,6)=6 

np(2,l)=2 

np(2,2)=7 

np(2,3)=8 

np(2,4)=9 

np(2,5)=10 

np(2,6)=ll 

np(3,l)=3 
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np(3,2)=8 

np(3,3)=12 

np(3,4)=13 

np(3,5)=14 

np(3,6)=15 

np(4,l)=4 

np(4,2)=9 

np(4,3)=13 

np(4,4)=16 

np(4,5)=17 

np(4,6)=18 

np(5,l)-5 

np(5,2)=10 

np(5,3)=14 

np(5,4)=17 

np(5,5)=19 

np(5,6)=20 

np(6,l)=6 

np(6,2)=ll 

np(6,3)=15 

np(6,4)=18 

np(6,5)=20 

np(6,6)=21 


1=0 

do  i=l,nofsps 

do  j=l,nofsps 
1=1+1 
mp(i,j)=l 
end  do 

end  do 

vmi(l)=2.0158 
vmi(2)=  17.03 17 
vmi(3)=l  14.8211 
vmi(4)=l  3 1.8528 
vmi(5)=347.4336 
vmi(6)=28.016 
c sigii(l)=2.968 
sigii(l)=2.915 
sigii(2)=3.4 
sigii(3)=5.47 
sigii(4)=5.89 
sigii(5)=7.631 
sigii(6)=3.652 
epii(l)=33.3 


c 
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epii(l)=38.0 

epii(2)=220. 

epii(3)=378.2 

epii(4)=278.9 

epii(5)=273.2 

epii(6)=98.4 

do  i=l,nofsps 

do  j=l,nofsps 

ij=np(i,j) 

sigij  (ij)=(sigii(i)+sigiiO  ))/2 . 
epkij  (ij  )=dsqrt(epii(i)  * epii(j)) 

end  do 

end  do 

sigij  (6)-3. 348 
epkij(6)=50.32 
sigij  (2)=3 .2093 
epkij  (2)=80.99 
sigij(ll)=3.5191 
epkij(l  1)=146.906 
do  j=l,non 

yitol=0. 

vmm(i)=0. 

c do  i=l,nofsps-l 

c yitol=yim+xl(nc(  1 05 0)+20- 1 +i+nn(j )) 

c end  do 

c if  (yitol  .ge.  1.0)  then 

c do  i=l,nofsps-l 

c yi(i,j)=xl(nc(1050)+20-l+i+nn(j))/yitol 

c normalize  yi 

c end  do 

c yi(nofsps,j)=0.0 

c assuming  last  is  main  carrier  gas 

c else 

c do  i=l,nofsps-l 

c yi(iJ)=xl(nc(1050)+20-l+i+nn(j)) 

c end  do 

c yi(nofsps,j)=1.0-yitol 

c end  if 

end  do 
do  i=l,non 

c doj=l,nofsps 

c yi(j,i)=l./nofsps 

c end  do 

yi(l,i)=8.9d-7 
yi(2,i)=.37626 


5^  R? 


yi(3,i)=5.073d-5 
yi(4,i)=.005826 
yi(5,i)=1.175d-5 
yi(6,i)=.6177 
end  do 
do  j=l,non 

do  i=l,nofsps 

vmm(j  )=vmm(j  )+yi(i,j  )/vmi(i) 

end  do 

vm(j)=l./vmm(j) 

end  do 
do  j=l,non 

do  i=l,nofsps 

xi(i,j)=yi(i,j)*vm(])/vmi(j) 

end  do 

end  do 

1 St  Step —assemble  collision  intergrals— — 


pi=3. 141 592654 
vk=1.380658d-23 
h2-n2 


npoint=6 
do  i=l,non 

omigall(i)=17.008048+(2.4667414d-6)*(dlog(at(i)))*at(i)- 
& 1 ,4422292*(dlog(at(i)))+1762.5095/at(i)/dsqrt(at(i)) 

omigal2(i)=15.747647-1.3435928*(dlog(at(i)))+4.0131919/ 

& dsqrt(at(i))+8 18.788  5 3/at(i)/dsqrt(at(i)) 

omigal  3(i)=l  7.786468+.037 1 02 1 44*((dlog(at(i)))*  *2)- 
1.927652*(dlog(at(i)))+ 
853.64569*(dlog(at(i)))/at(i)/at(i) 
omiga22(i)=- 16.1 42395+2 12.33759/(dlog(at(i)))-2 16.20479/ 

& dsqrt(at(i))+500.39094/at(i)/at(i) 

omijl  1(6, i)=omigal  l(i)/3. 348/3. 348 
omij  1 2(6,i)=omigal  2(i)/3 .348/3.348 
omijl  3(6, i)=omigal3(i)/3. 348/3. 348 
omiJ22(6,i)=omiga22(i)/3. 348/3. 348 


end  do 


h2-nh3 

h2-h2 


npoint=2 
sigi=sigii(l) 
do  i=l,non 

omigal  l(i)=14.941 186+.041 140446*((dlog(at(i)))**2)- 


Rp  Rp 
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1.7557228*(dlog(at(i)))+ 
916.39279*(dlog(at(i)))/at(i)/at(i) 
omigal2(i)=6.7374647+(5.1307243e-6)*at(i)*(dlog(at(i)))- 
& ,063643092*((dlog(at(i)))**2)+20.537264/dsqrt(at(i)) 

omigal3(i)=-11.404959+.0070803473*((dlog(at(i)))**2)+ 

& 123.191 68/(dlog(at(i)))-84.62 1 6 1 2/dsqrt(at(i)) 

omiga22(i)=14.006899-1.2119745*(dlog(at(i)))+55.40548/at(i)+ 

& 1381.1616/at(i)/at(i) 

omij  11(1  ,i)=omigal  1 (i)/sigi/sigi 
omij  1 2(1  ,i)=omigal  2(i)/sigi/sigi 
omij  13(1  ,i)=omiga  1 3 (i)/ sigi/sigi 
omiJ22(l,i)=omiga22(i)/sigi/sigi 
epi=epii(l) 
sigi=sigii(l) 

call  collint(epi,omigal  l,omigal2,omigal3,omiga22) 
omij  1 1 (1  ,i)=omigal  1 (i) 
omij  1 2(  1 ,i)=omigal  2(i) 
omij  13(1  ,i)=omiga  1 3 (i) 
omij  22(  1 ,i)=omiga22(i) 

end  do 

c nti3 

epi=epii(2) 

sigi=sigii(2) 

call  collint(epi,omigal  l,omigal2,omigal3,omiga22) 
c npoint=7 

do  i=l,non 

omij  1 1 (7,i)=omigal  1 (i) 
omij  1 2(7,i)=omiga  1 2(i) 
omij  1 3(7,i)=omigal  3(i) 
omij  22(7, i)=omiga22(i) 
c npoint=2 

omij  1 1 (2,i)=(omij  1 1 (1  ,i)+omij  1 1 (7,i))/2. 
omij  1 2(2,i)=(omij  1 2(  1 ,i)+omij  1 2(7,i))/2. 
omij  1 3(2,i)=(omij  13(1  ,i)+omij  1 3(7,i))/2. 
omij22(2,i)=(omij22(l,i)+omij22(7,i))/2. 

end  do 


c n2-nh3 

c n2-n2 

sigi=sigii(6) 
do  i=l,non 

omigal  1 (i)=9.9609368-.00 1 08 1 8748*at(i)+(  1 ,0806343e-7)*at(i)* 
& at(i)+749.1 1747/at(i)+6003.6332/at(i)/at(i) 

omigal  2(i)=4.8743454-.00441 67543  *(dlog(at(i)))*dsqrt(at(i))+ 

& 35.7661 58/(dlog(at(i)))+28621.337/at(i)/at(i) 
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omigal3(i)=9.7323604-.0014811822*at(i)+(1.8424564e-7)*at(i)* 
& at(i)+260.61415/at(i)+14646.726/at(i)/at(i) 

omiga22(i)=10.263407-.00042658832*at(i)-(2.2261337e-8)* 

& at(i)*at(i)+1094.0425/at(i)-4453.7957/at(i)/at(i) 

end  do 

c npoing=21 
do  i=l,non 

omijl  l(21,i)=omigal  l(i)/sigi/sigi 
omij  1 2(21  ,i)=omigal  2(i)/sigi/sigi 
omij  13(21  ,i)=omiga  1 3 (i)/sigi/ sigi 
omij  22(2 1 ,i)=omiga22(i)/sigi/ sigi 
epi=epii(6) 
sigi=sigii(6) 

call  collint(epi,omigal  l,omigal2,omigal3,omiga22) 
omij  11(21  ,i)=omiga  1 1 (i) 
omij  1 2(2 1 ,i)=omigal  2(i) 
omij  13(21  ,i)=omigal  3(i) 
omij22(21  ,i)=omiga22(i) 
c npoing=l  1 

omij  11(11  ,i)=(omij  1 1 (2 1 ,i)+omij  1 1 (7,i))/2. 
omij  12(11  ,i)=(omij  1 2(2 1 ,i)+omij  1 2(7,i))/2. 
omij  13(11  ,i)=(omij  13(21  ,i)+omij  1 3(7,i))/2. 
omij  22(11  ,i)=(oiTiij  22(2 1 ,i)+omij  22(7,i))/2 . 


■2nd  n2-h2 


omijl  1(6, i)=(omijl  1(21, i)+omijl  1(1, i))/2. 
omij  1 2(6,i)=(omij  1 2(2 1 ,i)+omij  1 2(  1 ,i))/2. 
omij  1 3(6,i)=(omij  1 3(2 1 ,i)+omij  1 3(1  ,i))/2. 
omij  22(6,i)=(omij  22(2 1 ,i)+omij22(  1 ,i))/2. 

end  do 
c 


c — h2-tmg — 

c tmg-tmg 

epi=epii(3) 

sigi=sigii(3) 

call  collint(epi,omigal  l,omigal2,omigal3,omiga22) 
c npoing=12  (tmg-tmg) 


do  i=l,non 

omij  11(1 2,i)=omigal  1 (i) 
omij  1 2(  1 2,i)=omigal  2(i) 
omij  13(1 2,i)=omiga  1 3 (i) 
omij  22(  1 2,i)=omiga22(i) 

end  do 

c npoing=3  (h2-tmg) 
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do  i=l,non 

omijl  1(3, i)==(ornijl  1(1, i)+omijl  1(12, i))/2. 
omij  1 2(3  ,i)=(omij  1 2(  1 ,i)+omij  1 2(  1 2,i))/2. 
omij  13(3,i)=(omij  1 3(l,i)+omij  1 3(12,i))/2. 
omij22(3  ,i)=(omij  22(  1 ,i)+omij  22(  1 2,i))/2. 

end  do 

^j|C**j|Cj|t***>l'>t!>t!!|!*******=|!!t!  + + ******!|‘=l==l==l==l=**********=l'****!|‘H!****H«>|!**H<**=l<!M<>HH<H=!M< 


c h2-tmg:nh3 

c adduct-adduct 

epi=epii(4) 

sigi=sigii(4) 

call  collint(epi,omigal  l,omigal2,omigal3,omiga22) 
c npoing- 16  (adduct-adduct) 


do  i=l,non 

omijl  l(16,i)=omigal  l(i) 
omij  1 2(1 6,i)=omigal  2(i) 
omij  13(1 6,i)=omigal  3 (i) 
omij22(  1 6,i)=omiga22(i) 

end  do 

c npoing=4  (h2-adduct) 
do  i=l,non 

omijl  1(4, i)=(omijl  1(1, i)+omijl  1(1 6,i))/2. 
omij  12(4,i)=(omij  12(1  ,i)+omij  12(1 6,i))/2. 
omij  1 3(4,i)=(omij  1 3(1  ,i)-^omij  1 3(1 6,i))/2. 
omij  22(4,i)=(omij  22(  1 ,i)+omij  22(1 6,i))/2 . 

end  do 


c h2-oligomer 

c oligomer-oligomer 

epi=epii(5) 

sigi=sigii(5) 

call  collint(epi,omigal  l,omigal2,omigal3,omiga22) 
c npoing=19  (oligomer-oligomer) 


do  i=l,non 

omijl  1(1 9,i)=omigal  l(i) 
omij  1 2(  1 9,i)=omigal  2(i) 
omij  1 3(1 9,i)=omigal  3(i) 
omij22(l  9,i)=omiga22(i) 

end  do 

c npoing=5  (h2-oligomer) 
do  i=l,non 

omijl  1(5, i)=(omijl  1(1, i)+omijl  1(1 9,i))/2. 
omij  1 2(5,i)=(omij  1 2(  1 ,i)+omij  1 2(  1 9,i))/2. 
omij  1 3(5,i)=(omij  1 3(1  ,i)+omij  1 3(19,i))/2. 
omij22(5,i)=(omij22(l,i)+omij22(19,i))/2. 
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end  do 

c 

**j|c*!|c**i|cJlc>|c:tci|=i|c**i|cHc**  + + ***>l'  + >l'****  + ***=M‘*=l=**************>l=****=l<=l=**=l<=l'=l'***=i=* 

L' 


c nh3-tmg- 

c npoing=8  (nh3-tmg) 


do  1=1, non 

omij  1 l(8,i)=(omij  1 l(7,i)+omij  1 l(12,i))/2. 
omij  12(8, i)=(omij  12(7, i)+omij  12(12, i))/2. 
omij  1 3(8,i)=(omij  1 3(7,i)+omij  1 30  2,i))/2. 
omij22(8,i)=(omij22(7,i)+omiJ22(12,i))/2. 

end  do 

c 


c nh3 -adduct 

c npoing=9  (nh3 -adduct) 


do  1=1, non 

omij  1 1 (9,i)=(omij  1 1 (7,i)-tomij  1 1 (1 6,i))/2. 
omij  1 2(9,i)=(omij  1 2(7,i)+omij  1 2(  1 6,i))/2. 
omij  1 3(9,i)=(omij  1 3(7,i)+omij  1 30  6,i))/2. 
omij22(9,i)=(omij22(7,i)-Homij22(16,i))/2. 

end  do 

c - - 


c nh3 -oligomer 

c npoing=10  (nh3 -oligomer) 


do  i=l,non 

omij  11(1 0,i)=(omij  1 1 (7,i)+omij  11(1 9,i))/2. 
omij  1 2(  1 0,i)=(omij  1 2(7,i)+omij  12(1 9,i))/2. 
omij  1 3(1 0,i)=(omij  1 3(7,i)+omij  13(1 9,i))/2. 
omij22(10,i)=(omij22(7,i)+omij22(19,i))/2. 

end  do 

c 

_*H«**H<**=t=H=>l<****=(=*!|=*************>l=>l<  + *****H”K*Hc*****>l<>l<>l<*=l==l==i!***=|c********>|s*** 


c tmg-adduct 

c npoing=13  (tmg-adduct) 


do  i=l,non 

omij  11(13  ,i)=(omij  11(1 2,i)-i-omij  11(1 6,i))/2 . 
omij  12(13  ,i)=(omij  1 2(  1 2,i)+omij  1 2(  1 6,i))/2. 
omij  1 3 (1 3 ,i)=(omij  1 3 (1 2 ,i)-tomij  1 3 ( 1 6,  i))/2 . 
omij22(13,i)=(omij22(12,i)+omij22(16,i))/2. 

end  do 

c 

_j|CjK!(!>|C  + **H'******=l‘*s|'*>|:!|!>(!****H‘************=l'*****>M‘>i'>l<>M'=l<>|!***=M'H'H'!|=s(!**»l<**H<!|'** 


c tmg-oligomer- 

c npoing=14  (tmg-oligomer) 
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do  i=l,non 

omijl  1(14, iMomijl  1(12, i)+omijl  1(1 9,i))/2. 
omij  1 2(  1 4,i)=(omij  1 2(  1 2,i)+omij  1 2(  1 9,i))/2. 
omij  13(1 4,i)=(omij  13(1 2,i)+omij  13(1 9,i))/2. 
omij22(14,i)==(omij22(12,i)+omij22(19,i))/2. 

end  do 

c — 


c - adduct-oligmer- 

c npoing=17  (adduct-oligomer) 


do  i=l,non 

omij  1 1(1 7,i)=(omijl  1(1 6,i)+omijl  1(1 9,i))/2. 
omij  12(1 7,i)=(omij  1 2(  1 6,i)+omij  1 2(  1 9,i))/2. 
omij  1 3(1 7,i)=(omij  1 3(1 6,i)+omij  13(1 9,i))/2. 
omij22(l  7,i)=(omij22(  1 6,i)+omij22(  1 9,i))/2. 

end  do 

c 


c - n2-tmg- 

c npoing=15  (n2-tmg) 


do  i=l,non 

omijl  l(15,i)=(omijl  l(21,i)+omijl  l(12,i))/2. 
omij  1 2(  1 5,i)=(omij  1 2(2 1 ,i)+omij  1 2(  1 2,i))/2. 
omij  1 3(1 5,i)=(omij  1 3(2 1 ,i)+omij  13(1 2,i))/2. 
omij22(15,i)=(omij22(21,i)+omij22(12,i))/2. 

end  do 

c 


c n2-adduct 

c npoing=18  (n2-adduct) 


do  i=l,non 

omij  11(1 8,i)=(omij  11(21  ,i)+omij  11(1 6,i))/2. 
omij  1 2(  1 8,i)=(omij  1 2(2 1 ,i)+omij  12(1 6,i))/2. 
omij  1 3(1 8,i)=(omij  13(21  ,i)+omij  13(1 6,i))/2. 
omij22(l  8,i)=(omij22(2 1 ,i)+omij22(  1 6,i))/2. 

end  do 

c 

^!|<!(!H<*HeJlCj|C)|e*  + *i|c  + **j|<>|S*H!****  + *=|!****=l==l'********=l‘=|!S(!**>|t>l‘*=l'>l«***=M«**s|!****!(!=M'***** 

c n2-oligomer 

c npoing=20  (n2-oligomer) 
do  i=l,non 

omijl  l(20,i)=(omijl  1(21, i)+omijl  1(1 9,i))/2. 
omij  1 2(20,i)=(omij  1 2(2 1 ,i)+omij  1 2(  1 9,i))/2. 
omij  1 3(20,i)=(omij  1 3(2 1 ,i)+omij  1 3(1 9,i))/2. 
omij 22(20,i)=(omij  22(2  l,i)+omij22(19,i))/2. 
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end  do 


c 

c diffusivity — 

call  multid 


C 

C 

c 

c 

c 

c 

c 

c 


& 


& 


return 

end 

################################################### 

##  SUBROUTINE  collint  ## 

##  FEM  VERSION  2.0  ## 

##  min  huang  ## 

################################################### 

SUBROUTINE  collint(epi,omigal  l,omigal2,omigal3,omiga22) 

###################################################  DECLARATION 
IMPLICIT  REAL* 8 (A-H,0-Z) 

common/integrl/omijl  1(55, 20), omij  12(55, 20), omij  13(55, 20), 

omij22(55,20) 

common/spsinx/nofsps,non,np(  1 0, 1 0),mp(  10,10) 
common/species/  vmi(10),yi(10,20),xi(10,20),epkij(55),sigij(55) 
common/dij  dti/dij  1 (5  5 ,20),dij  2(5  5 ,20),bdij  (5  5 ,20),dti(5  5 ,20) 
common/temp/at(20) 

dimension  omigal  I(20),omigal2(20),omigal3(20),omiga22(20) 
print  *,  "epi",epi 
do  i==l,non 

tstar=at(i)/epi 

omiga  1 1 (i)=dexp(.295402-. 5 1 0069 * dlog(tstar)+ 
.189395*((dlog(tstar))**2)- 
.045427*  ((dlog(tstar))*  * 3)+ 

.0037928*((dlog(tstar))*  *4)) 
omigal  2(i)=.829977+.  1 262633  *dIog(tstar)- 
.045427*((dlog(tstar))**2)+ 

.005057067*((dlog(tstar))**3) 
omiga  1 2(i)=omiga  1 2(i)  * omiga  1 1 (i) 
omigal  3(i)=.87 1 317631 +.0265709*(dlog(tstar))+ 
.010875045*((dlog(tstar))**3)- 
,018361235*((dlog(tstar))**4)+ 

.00583381  *((dlog(tstar))**5)- 
.000804046*((dlog(tstar))**6)+ 
.000044753*((dlog(tstar))**7) 
omiga  1 3 (i)=omigal  3 (i)  * omiga  1 2(i) 
omiga22(i)=dexp(.4664 1 -.5699 1 *dlog(tstar)+ 

.19591  *((dlog(tstar))**2)- 


n n 
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& ,03879*((dlog(tstar))**3)+ 

& .00259*((dlog(tstar))**4)) 

print  *,  omigal  l(i),omigal2(i),omigal3(i),omiga22(i) 
end  do 
return 
end 
C 

c ################################################### 

C ##  SUBROUTINE  multid  ## 

C ##  FEM  VERSION  2.0  ## 

C ##  min  huang  ## 

C ################################################### 

c 

SUBROUTINE  multid 


###################################################  DECLARATION 
IMPLICIT  REAL* 8 (A-H,0-Z) 

common/integrl/omij  I I(55,20),omij  12(55,20),omij  1 3(55,20), 

& omij22(55,20) 

common/spsinx/nofsps,non,np(  1 0, 1 0),mp(  1 0, 1 0) 

common/species/  vmi(I0),yi(10,20),xi(10,20),epkij(55),sigij(55) 

common/dijdti/dijl(55,20),dij2(55,20),bdij(55,20),dti(55,20) 

common/temp/at(20) 

common/pred/cij  (5  5 ,20) 

dimension  vki2(5  5 ,20),astar(  10,1 0,20),bstar(  1 0, 1 0,20) 

& ,cstar(  1 0, 1 0,20),slamij(  1 0, 1 0,20),clam00(  10,10,20) 

& ,clam01(10,10,20),claml  1(1 0,1 0,20),detij2(l  0,1 0,20) 

dimension  clam  10(10,1 0,20),vmi2(  1 0),pdet  1(10,1 0),pdet2(  10,10) 

& ,cl  1 1 ( 1 0, 1 0),cliv(  1 0, 1 0),cliv  11(10,1 0,20) 

& ,cl0 1(10,1 0),cl  1 0(  1 0, 1 0),zero(  1 0, 1 0),vki  1(55 ,20) 

print  *,  "4" 
pi=3. 141 592654 
vk=1.380658d-23 
avn=6.022d23 
do  i-l,nofsps 
vmi2(i)=  vmi(i)/avn/1000. 
end  do 


write  (1,*)  "k,j,slamij(k,j,l)=" 
do  i=l,non 
c kj=0 

do  k=l,nofsps 

do  j=l,nofsps 
kj=np(k,j) 

astar(k,j  ,i)=omij  22(kj  ,i)/omij  1 1 (kj  ,i) 
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bstar(k,j,i)=(5.*omij  1 2(kj  ,i)-4.  *omij  1 3(kj,i))/omij  1 1 (kj  ,i) 
cstar(k,j  ,i)=omij  1 2(kj  ,i)/omij  1 1 (kj  ,i) 

c common  const,  parameter  was  removed  for  reduce  error- 

slamij  (k,j  ,i)=  1.9891  d-2*  4. 1 84  * (dsqrt((vmi(k)+vmi(j )) 

& /vmi(k)/vmi(j)/2.*at(i))) 

& /omij22(kj,i)/sigij(kj)**2 

c slamij  (k,j,i)=  75./32.*(dsqrt((vmi2(k)+vmi2(j)) 
c & /vmi2(k)/vmi2(j)/8./pi*at(i)*vk**3)) 

c & /omij22(kj,i)/(l.d-10*sigij(kj))**2 

write  (1,*)  k,j,slamij(k,j,l) 
end  do 

end  do 


end  do 
do  i=l,non 

do  k=l,nofsps 

do  j=l,nofsps 

clam00(k,j,i)=0-0 

clam01(kj,i)=0.0 
claml  l(k,j,i)=0-0 

end  do 

end  do 

end  do 
do  i=l,non 

do  k=l,nofsps 


do  j=l,nofsps 

kj2=mp(k,j) 
kj=np(k,j) 
if  (j  .ne.  k)  then 
clam00(k,j  ,i)=-xi(k,i)*xi(j  ,i) 

& /2./astar(k,j,i)/slamij(k,j,i) 

clam01(k,j,i)=+xi(k,i)*xi(j,i)*(6.*cstar(k,j,i)-5.)* 

& vmi2(k)/(vmi2(k)+vmi2(j))/astar(k,j,i)/slamij(k,j,i)/4. 

claml  I(k,j,i)=+xi(k,i)*xi(j,i)*vmi2(k)*vmi2(j)* 

& (4 . * astar(k,j  ,i)+3 . * bstar(k,j  ,i)-5  5 ./4.)/astar(k,j  ,i) 

& /slamij  (k,j  ,i)/(vmi2(k)+vmi2(j  ))/(vmi2(k)+vmi2(j  ))/2 . 

c 

C k=j 

else 

clam00(k,j,i)=0.0 

clam01(k,j,i)=0-0 

claml  l(k,j,i)=+xi(k,i)*xi(k,i)/slamij(k,j,i) 

c summation  part — — 

do  l=l,nofsps 

ll=np(k,l) 
if  (1  .ne.  k)  then 


Rp  R?  Rp  Rp  Rp  Rp 
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claml  l(k,j,i)=claml  l(k,j,i)+xi(k,i)*xi(l,i) 

*(7.5*vmi2(k)*vmi2(k)+6.25*vmi2(l)*vmi2(l)-3. 

*vmi2(l)*vmi2(l)*bstar(k,l,i)+4.*vmi2(k)*vmi2(l) 

*astar(k,l,i))/(vmi2(k)+vmi2(l))/(vmi2(k) 

+vmi2(l))/astar(k,l,iyslamij(k,l,iy2. 


clamO  1 (k,j  ,i)=clamO  1 (k,j,i)-xi(k,i)*xi(l,i) 

*(6.*cstar(k,l,i)-5.)*vmi2(l) 

/(vmi2(k)+vmi2(l)yastar(k,l,iyslamij(k,I,i)/4. 


clamOO(k,j,i)=clamOO(k,j,i)+xi(k,i)*xi(l,i) 

& /astar(k,l,i)/slamij(k,l,iy2. 

end  if 

30  end  do 

end  if 

end  do 

end  do 

do  k=l,nofsps 

do  j=l,nofsps 

kj2=mp(k,j) 
if  (j  .ne.  k)  then 

clam  1 0(k,j  ,i)=clam0 1 (k,j  ,i)*  vmi2(j)/vmi2(k) 

else 

clam  1 0(k,j  ,i)=clam0 1 (k,j  ,i) 

end  if 

end  do 

end  do 

do  k=l,nofsps 

do  j=l,nofsps 

kj2=mp(k,j) 
if  (i  .ne.  k)  then 

clam  1 0(k,j  ,i)=clam0 1 (k,j  ,i)  * vmi(j )/ vmi(k) 

else 

clam  1 0(k,j  ,i)=clam0 1 (k,j  ,i) 

end  if 

end  do 

end  do 

do  k=l,nofsps 

do  j=l,nofsps 

kj2=mp(k,j) 

write  ( 1 ,*)  j ,clam  1 0(k,j ,i),clam0 1 (k,j ,i) 

end  do 

end  do 

do  k=l,nofsps 

do  j=l,nofsps 


147 


kj2=mp(kj) 

clll(k,j)=clamll(k,j,i) 

cliv(k,j)=0.0 

zero(k,j)=0.0 

write  (1,*)  "clamll(k,j,i)=",k,j,clamll(kj,i) 
end  do 

end  do 

call  MATOPR2  ('INV',cll  l,zero,cliv,nofsps) 

do  k=l,nofsps 

do  j=l,nofsps 

kj2=mp(k,j) 
clivl  l(k,j,i)=cliv(k,j) 

write  (1,*)  "clivl  l(k,j,i)=",k,j,clivl  l(kj,i) 
end  do 

end  do 

end  do 

do  i=l,nofsps 
cl  0=0. 

do  j=l,nofsps 

c 1 0=c  1 0+clam  1 0(i,j , 1 ) 

end  do 

write  (1,*)  "cl0=",cl0 

end  do 

do  j=l,nofsps 
c01=0. 

do  i=l,nofsps 

cO  1 =c0 1 +clam0 1 (ij , 1 ) 

end  do 

write  (1,*)  "c01=",c01 

end  do 

c deltaij  2nd  - 

do  i=l,non 

do  k=l,nofsps 

do  j=l,nofsps 

kj2=mp(k,j) 

detij2(k,j,i)=.0 

end  do 

end  do 

do  k=l,nofsps 

do  j=l,nofsps 

kj2=mp(k,j) 
cliv(kj)=clivl  l(k,j,i) 
clO  1 (k,j)=clam0 1 (k,j  ,i) 
cl  1 0(k,j)=clam  1 0(k,j  ,i) 

end  do 
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end  do 

call  MATOPR2  ('MUL’,cliv,cllO,pdetl,nofsps) 
call  MATOPR2  ('MUL',cl01,pdetl,pdet2,nofsps) 
do  k=l,nofsps 

do  j=l,nofsps 

detij  2(k,j  ,i)=pdet2(k,j  )/clam00(k,j  ,i) 
write  (1,*)  "detij2(k,j,i)=",kj,detij2(k,j,i) 
end  do 

end  do 

end  do 

c binary  diffusivity  1 ,2 

do  i=l,non 

do  k=l,nofsps 

do  j=l,nofsps 

kj-np(kj) 
if  (i  .ne.  k)  then 

vv=(vmi(k)+vnii(j  ))/2  ./vmi(k)/vmi(j  )*  ((at(i))*  * 3 ) 
dijl(kj,i)=.002628*dsqrt(w)/omijl  l(kj,i)/sigij(kj)/sigij(kj)/10000. 
fs=detij2(kj,i) 
dij2(kj,i)=(dij  l(kj,i))/(l  .-fs) 
else 

dijl(kj,i)=0.0 

dij2(kj,i)=0.0 

end  if 

write  (1,*)  dijl(kj,i),detij2(k,j,i),dij2(kj,i) 
end  do 

end  do 

end  do 

c 

c multi  diffusivity 

do  i=l,non 

kj=0 

do  k=l,nofsps 

do  j=k,nofsps 
kj-kj+1 

if  (j  .ne.  k)  then 

cij  (kj  ,i)=xi(k,i)*  xiO  ,i)/dij  1 (kj  ,i) 
else 

cij(kj,i)=.0 
end  if 


50  end  do 

end  do 

end  do 
call  obbdij 

c multi  thermal  diff- 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


& 


do  i=l,non 

do  k=l,nofsps 

do  j=l,nofsps 

kj2=mp(k,j) 

cll  l(k,j)=claml  l(k,j,i) 

cliv(k,j)=0.0 

zero(k,j)=0.0 

end  do 

end  do 

call  MATOPR2  ('INV',cll  l,zero,cliv,nofsps) 
do  k=l,nofsps 

do  j=l,nofsps 

kj2=mp(k,j) 
clivl  l(k,j,i)=cliv(k,j) 

end  do 

end  do 

do  k=l,nofsps 

vki2(k,i)=.0 
dti(k,i)=.0 
do  l=l,nofsps 

vkil(l,i)=.0 
do  j=l,nofsps 

vki  1 (l,i)=vki  1 (l,i)+2.5  *xi(j  ,i)*cliv  1 1 (j  ,l,i) 

end  do 

end  do 

do  l=l,nofsps 

vki2(k,i)=vki2(k,i)+vki  1 (l,i)*  clam  1 0(1, k,i) 

end  do 

do  j=l,nofsps 

do  l=l,nofsps 

vki2(k,i)=vki2(k,i)+2.5*xi(j,i)*clam01(k,l,i) 

*clivll(l,j,i) 

end  do 

end  do 

end  do 

end  do 
vk=0. 

do  i=l,nofsps 

vk=vk+vki2(i,l) 

print  *,  "vki2(k,i)",vki2(i,l) 

write  (1,*)  "vki2(k,i)",vki2(i,l) 

end  do 

print  *,  "sig-vki=",vk 
write  (1,*)  "sig  vki=",vk 


n o 
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c 

c 


C 

C 

C 

C 

c 

c 

c 


DTi 

do  i=l,non 

do  k=l,nofsps 

do  j=l,nofsps 

l=np(k,j) 

dti(k,i)=dti(k,i)+bdij(l,i)*vki2(i,i) 

end  do 

end  do 

end  do 

alfa 

return 

end 

################################################### 

##  SUBROUTINE  obbdij  ## 

##  FEM  VERSION  2.0  ## 

##  min  huang  ## 

################################################### 

SUBROUTINE  obbdij 

###################################################  DECLARATION 

IMPLICIT  REAL* 8 (A-H,0-Z) 

common/spsinx/nofsps,non,np(  1 0, 1 0),mp(  10,10) 

common/species/  vmi(  1 0),yi(  1 0,20), xi(  1 0,20),epkiJ  (55),sigij  (55) 

common/dijdti/dijl(55,20),diJ2(55,20),bdij(55,20),dti(55,20) 

commonytemp/at(20) 

common/pred/cij  (5  5 ,20) 

dimension  c(55),pa(55,55),pai(55,55),pf(15),pd(15),bdii(10) 

& ,res(55,55) 

do  i=l,non 
kj=0 

do  k=l,nofsps 

do  j=k,nofsps 
kj=kj+l 
c(kj)=cij(kj,i) 
bdij(kj,i)=0. 

end  do 

end  do 

dokl=l,15 

pf(kl)=0. 
do  k2=l,15 
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pa(kl,k2)=0. 

end  do 

end  do 


c build  parameter  matrix 

pa(0 1 ,0 1 )=-yi(  1 ,i)*(2.  *C(02)+c(08)+c(09)+c(  1 0)+c(  1 1 )) 
& -yi(2,i)*(2.*c(02)+c(03)+c(04)+c(05)+c(06)) 

pa(01,02)=+yi(l,i)*c(08)-yi(3,i)*c(2) 
pa(01,03)=+yi(l,i)*c(09)-yi(4,i)*c(2) 
pa(0 1 ,04)=+yi(l  ,i)*c(l  0)-yi(5,i)*c(2) 
pa(0 1 ,05)=+yi(  1 ,i)*c(  1 1 )-yi(6,i)*c(2) 
pa(01,06)=+yi(2,i)*c(03)-yi(3,i)*c(2) 
pa(01,07)=+yi(2,i)*c(04)-yi(4,i)*c(2) 
pa(01,08)=+yi(2,i)*c(05)-yi(5,i)*c(2) 
pa(01,09)=+yi(2,i)*c(06)-yi(6,i)*c(2) 
pf(l)  =2.*yi(l,i)*yi(2,i) 


pa(02,02)=-yi(l,i)*(2.*C(03)+c(08)+c(13)+c(14)+c(15)) 
& -yi(3,i)*(c(02)+2.*c(03)+c(04)+c(05)+c(06)) 

pa(02,01)=+yi(l,i)*c(08)-yi(2,i)*c(3) 
pa(02,03)=+yi(l,i)*c(13)-yi(4,i)*c(3) 
pa(02,04)=+yi(l,i)*c(14)-yi(5,i)*c(3) 
pa(02,05)=+yi(l,i)*c(15)-yi(6,i)*cP) 
pa(02,06)=+yi(3,i)*c(02)-yi(2,i)*c(3) 
pa(02,10)=+yi(3,i)*c(04)-yi(4,i)*c(3) 
pa(02, 1 1 )=+yi(3  ,i)*  c(05)-yi(5  ,i)  * c(3 ) 
pa(02,12)=+yi(3,i)*c(06)-yi(6,i)*c(3) 
pf(2)  =2.*yi(l,i)*yi(3,i) 


pa(03 ,03)=-yi(  1 ,i)*  (2.  *C(04)+c(09)+c(  1 3)+c(  1 7)+c(  1 8)) 
& -yi(4,i)*(C(02)+c(03)+2.*c(04)+c(05)+c(06)) 

pa(03 ,0 1 )=+yi(  1 ,i)*  c(09)-yi(2,i)*  c(4) 
pa(03,02)=+yiO  ,i)*c(  1 3)-yi(3,i)*c(4) 
pa(03,04)=+yi(  1 ,i)*c(  1 7)-yi(5,i)*  c(4) 
pa(03,05)=+yi(l,i)*c(18)-yi(6,i)*c(4) 
pa(03,07)=+yi(4,i)*c(02)-yi(2,i)*c(4) 
pa(03,10)=+yi(4,i)*c(03)-yi(3,i)*c(4) 
pa(03, 1 3)=+yi(4,i)*c(05)-yi(5,i)*c(4) 
pa(03,14)=+yi(4,i)*c(06)-yi(6,i)*c(4) 
pf(3)  =2.*yi(l,i)*yi(4,i) 


pa(04,04)=-yi(  1 ,i)*  (2.  * C(05)+c(  1 0)+c(  1 4)+c(  1 7)+c(20)) 
& -yi(5,i)*(c(02)+c(03)+c(04)+2.*c(05)+c(06)) 

pa(04,0 1 )=+yi(  1 ,i)*c(  1 0)-yi(2,i)*c(5) 
pa(04,02)=+yi(  1 ,i)*c(  1 4)-yi(3  ,i)*  c(5) 
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pa(04,03)=+yi(  1 ,i)*c(  1 7)-yi(4,i)*c(5) 
pa(04,05)=+yi(l,i)*c(20)-yi(6,i)*c(5) 
pa(04,08)=+yi(5,i)*c(02)-yi(2,i)*c(5) 
pa(04,ll)=+yi(5,i)*c(03)-yi(3,i)*c(5) 
pa(04, 1 3)=+yi(5,i)*c(04)-yi(4,i)*c(5) 
pa(04,15)=+yi(5,i)*c(06)-yi(6,i)*c(5) 
pf(4)  =2.*yi(l,i)*yi(5,i) 


pa(05,05)=-yi(  1 ,i)*(2.  *C(06)+c(  1 1 )+c(  1 5)+c(  1 8)+c(20)) 
& -yi(6,i)*  (c(02)+c(03)+c(04)+c(05)+2.  * c(06)) 

pa(05 ,0 1 )=+yi(  1 ,i)*c(  1 1 )-yi(2,i)*c(6) 
pa(05,02)=+yi(l,i)*c(15)-yi(3,i)*c(6) 
pa(05,03)=+yi(l,i)*c(18)-yi(4,i)*c(6) 
pa(05,04)=+yi(l,i)*c(20)-yi(5,i)*c(6) 
pa(05,09)=+yi(6,i)*c(02)-yi(2,i)*c(6) 
pa(05,12)=+yi(6,i)*c(03)-yi(3,i)*c(6) 
pa(05,14)=+yi(6,i)*c(04)-yi(4,i)*c(6) 
pa(05, 1 5)=+yi(6,i)*c(05)-yi(5,i)*c(6) 
pf(5)  =2.*yi(l,i)*yi(6,i) 


pa(06,06)=-yi(2,i)*(C(03)+2.*c(08)+c(13)+c(14)+c(15)) 
& -yi(3,i)*(c(02)+2.*c(08)+c(09)+c(10)+c(l  1)) 

pa(06,0 1 )=+yi(2,i)*c(03)-yi(  1 ,i)*  c(8) 
pa(06,07)=+yi(2,i)*c(13)-yi(4,i)*c(8) 
pa(06,08)=+yi(2,i)*c(14)-yi(5,i)*c(8) 
pa(06,09)=+yi(2,i)*c(15)-yi(6,i)*c(8) 
pa(06,02)=+yi(3  ,i)*c(02)-yi(  1 ,i)*c(8) 
pa(06, 1 0)=+yi(3  ,i)*c(09)-yi(4,i)*  c(8) 
pa(06,ll)=+yi(3,i)*c(10)-yi(5,i)*c(8) 
pa(06, 1 2)=+yi(3,i)*c(  1 1 )-yi(6,i)*c(8) 
pf(6)  =2.*yi(2,i)*yi(3,i) 


pa(07,07)=-yi(2,i)*(C(04)+2.*c(09)+c(13)+c(17)+c(18)) 

& -yi(4,i)*(c(02)+c(08)+2.*c(09)+c(10)+c(l  1)) 

pa(07,0 1 )=+yi(2,i)*c(04)-yi(  1 ,i)*c(9) 
pa(07,06)=+yi(2,i)*c(13)-yi(3,i)*c(9) 
pa(07,08)=+yi(2,i)*c(17)-yi(5,i)*c(9) 
pa(07,09)=+yi(2,i)*c(18)-yi(6,i)*c(9) 
pa(07,03)=+yi(4,i)*c(02)-yi(l,i)*c(9) 
pa(07,10)=+yi(4,i)*c(08)-yi(3,i)*c(9) 
pa(07,13)=+yi(4,i)*c(10)-yi(5,i)*c(9) 
pa(07,14)=+yi(4,i)*c(ll)-yi(6,i)*c(9) 
pf(7)  -2.*yi(2,i)*yi(4,i) 


pa(08,08)=-yi(2,i)*(C(05)+2.*c(10)+c(14)+c(17)+c(20)) 


-yi(5,i)*(c(02)+c(08)+c(09)+2.  *c(  1 0)+c(l  1 )) 
pa(08,01)=+yi(2,i)*c(05)-yi(l,i)*c(10) 
pa(08,06)=+yi(2,i)*c(14)-yi(3,i)*c(10) 
pa(08,07)=+yi(2,i)*c(17)-yi(4,i)*c(10) 
pa(08,09)=+yi(2,i)*c(20)-yi(6,i)*c(10) 
pa(08,04)=+yi(5,i)*c(02)-yi(l,i)*c(10) 
pa(08,l  l)=+yi(5,i)*c(08)-yi(3,i)*c(10) 
pa(08 , 1 3 )=+yi(5  ,i)*  c(09)-yi(4,i)*  c(  1 0) 
pa(08, 1 5)=+yi(5,i)*c(  1 1 )-yi(6,i)*c(  1 0) 
pf(8)  =2.*yi(2,i)*yi(5,i) 


pa(09,09)=-yi(2,i)*  (C(06)+2.  *c(  1 1 )+c(  1 5)+c(  1 8)+c(20)) 
-yi(6,i)*  (c(02)+c(08)+c(09)+c(  1 0)+2 . * c(  1 1 )) 
pa(09,0 1 )=+yi(2,i)*c(06)-yi(l  ,i)*c(l  1 ) 
pa(09,06)=+yi(2,i)*c(15)-yi(3,i)*c(l  1) 
pa(09,07)=+yi(2,i)*c(l  8)-yi(4,i)*c(  1 1 ) 
pa(09,08)=+yi(2,i)*c(20)-yi(5,i)*c(l  1) 
pa(09,05)=+yi(6,i)*c(02)-yi(l,i)*c(l  1) 
pa(09, 1 2)-+yi(6,i)*c(08)-yi(3,i)*c(  1 1 ) 
pa(09,14)=+yi(6,i)*c(09)-yi(4,i)*c(ll) 
pa(09,15)=+yi(6,i)*c(10)-yi(5,i)*c(l  1) 
pf(9)  =2.*yi(2,i)*yi(6,i) 


pa(10,10)=-yi(3,i)*(C(04)+c(09)+2.*c(13)+c(17)+c(18)) 
-yi(4,i)*  (c(03)+c(08)+2.  *c(  1 3)+c(  1 4)+c(  1 5)) 
pa(10,02)=+yi(3,i)*c(04)-yi(l,i)*c(13) 
pa(10,06)=+yi(3,i)*c(09)-yi(2,i)*c(13) 
pa(  10,11  )=+yi(3  ,i)*c(  1 7)-yi(5,i)*c(  1 3) 
pa(l  0, 1 2)=+yi(3,i)*c(l  8)-yi(6,i)*c(  13) 
pa(  1 0,03)=+yi(4,i)*c(03)-yi(l  ,i)*c(  13) 
pa(  1 0,07)=+yi(4,i)*c(08)-yi(2,i)*  c(  1 3) 
pa(  1 0, 1 3)=+yi(4,i)*c(  1 4)-yi(5  ,i)*  c(  1 3) 
pa(  1 0, 1 4)=+yi(4,i)*c(  1 5)-yi(6,i)*  c(  1 3) 
pf(10)  =2.*yi(3,i)*yi(4,i) 


pa(  1 1 , 1 1 )=-yi(3  ,i)*  (C(05  )+c(  1 0)+2 . * c(  1 4)+c(  1 7)+c(20)) 
-yi(5,i)*(c(03)+c(08)+c(  1 3)+2.  *c(  1 4)+c(  1 5)) 
pa(  1 1 ,02)=+yi(3,i)*c(05)-yi(  1 ,i)*c(  1 4) 
pa(  1 1 ,06)=+yi(3  ,i)*c(  1 0)-yi(2,i)*  c(  1 4) 
pa(  1 1 , 1 0)=+yi(3,i)*c(  1 7)-yi(4,i)*c(  1 4) 
pa(  1 1 , 1 2)=+yi(3  ,i)*c(20)-yi(6,i)  * c(  1 4) 
pa(  1 1 ,04)=+yi(5  ,i)  * c(03  )-yi(  1 ,i)  * c(  1 4) 
pa(  1 1 ,08)=+yi(5  ,i)*c(08)-yi(2,i)*  c(  1 4) 
pa(  1 1 , 1 3)==+yi(5  ,i)  *c(  1 3)-yi(4,i)  * c(  1 4) 
pa(ll,15)=+yi(5,i)*c(15)-yi(6,i)*c(14) 
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pf(ll)  =2.*yi(3,i)*yi(5,i) 


pa(  1 2, 1 2)=-yi(3,i)*(C(06)+c(  1 1 )+2.  *c(  1 5)+c(  1 8)+c(20)) 
& -yi(6,i)*(c(03)+c(08)+c(  1 3)+c(  1 4)+2.  * c(  1 5)) 

pa(  1 2,02)=+yi(3  ,i)*c(06)-yi(  1 ,i)*  c(  1 5) 
pa(12,06)=+yi(3,i)*c(ll)-yi(2,i)*c(15) 
pa(12,10)=+yi(3,i)*c(18)-yi(4,i)*c(15) 
pa(l  2, 1 1 )=+yi(34)*c(20)-yi(5,i)*c(  1 5) 
pa(12,05)=+yi(6,i)*c(03)-yi(l,i)*c(15) 
pa(  1 2,09)=+yi(6,i)*c(08)-yi(2,i)*c(  1 5) 
pa(  1 2, 1 4)=+yi(6,i)*c(  1 3)-yi(4,i)*  c(  1 5) 
pa(  1 2, 1 5)=+yi(6,i)*c(  1 4)-yi(5  ,i)*  c(  1 5) 
pf(12)  =2.*yi(3,i)*yi(6,i) 


pa(  1 3 , 1 3)=-yi(4,i)*(C(05)+c(  1 0)+c(  1 4)+2.  *c(  1 7)+c(20)) 
-yi(5,i)*(c(04)+c(09)+c(13)+2.*c(17)+c(18)) 
pa(  1 3 ,03)=+yi(4,i)*c(05)-yi(  1 ,i)*  c(  1 7) 
pa(  1 3 ,07)=+yi(4,i)*  c(  1 0)-yi(2,i)  * c(  1 7) 
pa(13,10)=+yi(4,i)*c(14)-yi(3,i)*c(17) 
pa(  1 3 , 1 4)=+yi(4,i)*c(20)-yi(6,i)*c(  1 7) 
pa(  1 3,04)=+yi(5,i)*c(04)-yi(l  ,i)*c(  1 7) 
pa(  1 3,08)=+yi(5,i)*c(09)-yi(2,i)*c(  1 7) 
pa(  1 3 , 1 1 )=+yi(5  ,i)  *c(  1 3 )-yi(3  ,i)  * c(  1 7) 
pa(  1 3 , 1 5)=+yi(5  ,i)  * c(  1 8)-yi(6,i)  * c(  1 7) 
pf(13)  =2.*yi(4,i)*yi(5,i) 


pa(l  4, 1 4)=-yi(4,i)*(C(06)+c(  1 1 )+c(  1 5)+2.*c(  1 8)+c(20)) 
& -yi(6,i)*(c(04)+c(09)+c(13)+c(17)+2.*c(18)) 

pa(14,03)=+yi(4,i)*c(06)-yi(l,i)*c(18) 
pa(  1 4,07)=+yi(4,i)*c(  1 1 )-yi(2,i)*  c(  1 8) 
pa(l  4, 1 0)=+yi(4,i)*c(l  5)-yi(3,i)*c(  1 8) 
pa(14,13)=+yi(4,i)*c(20)-yi(5,i)*c(18) 
pa(14,05)=+yi(6,i)*c(04)-yi(l,i)*c(18) 
pa(l  4,09)=+yi(6,i)*c(09)-yi(2,i)*c(  1 8) 
pa(14,12)=+yi(6,i)*c(13)-yi(3,i)*c(18) 
pa(  1 4, 1 5)=+yi(6,i)*c(  1 7)-yi(5,i)*c(  1 8) 
pf(14)  =2.*yi(4,i)*yi(6,i) 


pa(15,15)=-yi(5,i)*(C(06)+c(H)+c(15)+c(18)+2.*c(20)) 
-yi(6,i)*  (c(05)+c(  1 0)+c(  1 4)+c(  1 7)+2.  *c(20)) 
pa(l  5,04)=+yi(5,i)*c(06)-yi(l  ,i)*c(20) 
pa(15,08)=+yi(5,i)*c(l  l)-yi(2,i)*c(20) 
pa(  1 5, 1 1 )=+yi(5,i)*c(l  5)-yi(3,i)*c(20) 
pa(15,13)=+yi(5,i)*c(18)-yi(4,i)*c(20) 
pa(l  5,05)=+yi(6,i)*c(05)-yi(l  ,i)*c(20) 


& 
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pa(15,09)=+yi(6,i)*c(10)-yi(2,i)*c(20) 
pa(  1 5 , 1 2)=+yi(6,i)*  c(  1 4)-yi(3  ,i)*  c(20) 
pa(  1 5, 1 4)=+yi(6,i)*c(l  7)-yi(4,i)*c(20) 
pf(15)  =2.*yi(5,i)*yi(6,i) 

c 

c invert  pa 

call  gaussj  (pa,  15,15  ,pf,  1,1) 
do  n=l,15 
write  (1,*)  n,pf(n) 
end  do 

kj=0 

1=0 

do  k=l,nofsps 

bdii(k)=0.0 
do  j=k,nofsps 
kj=kj+l 
bdij(kj,i)=0.0 
if  (k  .eq.  j)  goto  10 
1=1+1 

bdij(kj,i)=pf(l) 


1 0 end  do 

end  do 

end  do 

c 

c find  Dii  from  eq.  1 0 pavlov 


do  i=l,non 

do  k=l,nofsps 

bdii(k)=0.0 
do  j=l,nofsps 

l=np(k,j) 

if  (j  .ne.  k)  bdii(k)=bdii(k)+yi(j,i)*bdij(l,i)/yi(k,i) 

end  do 

end  do 
kj=0 
1=0 

do  k=l,nofsps 

do  j=k,nofsps 
kj=kj+l 

if  (k  .eq.  j)  bdij(kj,i)=-bdii(k) 
end  do 

end  do 

end  do 

return 

end 

C ################################################### 


o n 
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C ##  SUBROUTINE  gaussj  ## 

C ##  FEM  VERSION  2.0  ## 

C ##  min  huang  ## 

C ################################################### 

c 

SUBROUTINE  gaussj  (a, n,np,b,m,mp) 

###################################################  DECLARATION 
IMPLICIT  REAL* 8 (A-H,0-Z) 
dimension  indxc(55),indxr(55),ipiv(55),a(55,55),b(l,55) 
do  j=l,n 

ipiv(i)=0 

enddo 
do  i=l,n 

big=0. 
do  j=l,n 

if(ipiv(j)  .ne.  1)  then 
do  k=l,n 

if  (ipiv(k)  .eq.  0)  then 

if  (dabs(a(j,k))  .ge.  big)  then 
big=dabs(a(j,k)) 
irow=j 
icol=k 

end  if 

else  if  (ipiv(k)  .gt.  1)  then 

pause  'sigular  matrix  in  gaussj' 

end  if 

enddo 

end  if 

enddo 

ipiv(icol)=ipiv(icol)+ 1 
if  (irow  .ne.  icol)  then 
do  1=1, n 

dum=a(irow,l) 
a(irow,l)=a(icol,l) 
a(icol,l)=dum 

enddo 
do  1=1, m 

dum=b(irow,l) 
b(irow,l)=b(icol,l) 
b(icol,l)=dum 

enddo 

endif 

indxr(i)=irow 
indxc(i)=icol 
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if  (a(icol,icol)  .eq.  0.)  pause  'singular  matrix  in  gaussj' 
pivinv=l  ./a(icol,icol) 
a(icol,icol)=l. 
do  1=1, n 

a(icol,l)=a(icol,l)*pivinv 

enddo 
do  1=1, m 

b(icol,l)=b(icol,l)*pivinv 

end  do 
do  ll=l,n 

if  (11  .ne.  icol)  then 

dum=a(ll,icol) 
a(ll,icol)=0. 
do  1=1, n 

a(ll,l)=a(ll,l)-a(icol,l)  * dum 

enddo 
do  1=1, m 

b(ll,l)=b(ll,l)-b(icol,l)*dum 

enddo 

endif 

end  do 

end  do 
do  l=n,l,-l 

if(indxr(l)  .ne.  indxc(l))  then 
do  k=l,n 

dum=a(k,indxr(l)) 

a(k,indxr(l))=a(k,indxc(l)) 

a(k,indxc(l))=dum 

enddo 

endif 

enddo 

return 

end 
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